# Inviscid Water-Waves and interface modeling

Emmanuel Dormy, Christophe Lacave

► To cite this version:

Emmanuel Dormy, Christophe Lacave. Inviscid Water-Waves and interface modeling. 2023. hal-04119856v3

HAL Id: hal-04119856

https://hal.science/hal-04119856v3

Preprint submitted on 12 Dec 2023

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
INVISCID WATER-WAVES AND INTERFACE MODELING

EMMANUEL DORMY AND CHRISTOPHE LACAVE

ABSTRACT. We present a rigorous mathematical analysis of the modeling of inviscid water waves. The free-surface is described as a parametrized curve. We introduce a numerically stable algorithm which accounts for its evolution with time. The method is shown to converge using approximate solutions, such as Stokes waves and Green-Naghdi solitary waves. It is finally tested on a wave breaking problem, for which an odd-even coupling suffices to achieve numerical convergence up to the splash without the need for additional filtering.

CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>1.1. Problem formulation</td><td>2</td></tr><tr><td>1.2. Numerical strategies</td><td>3</td></tr><tr><td>1.3. Objectives of the present study</td><td>4</td></tr><tr><td>1.4. Plan of the paper</td><td>5</td></tr><tr><td>2. Singular integrals representation at fixed time</td><td>6</td></tr><tr><td>2.1. Stream and potential functions</td><td>6</td></tr><tr><td>2.2. Potential and dipole formulae</td><td>10</td></tr><tr><td>2.3. Discussion about u<sub>ω,γ</sub> for the dipole formulation.</td><td>18</td></tr><tr><td>3. Evolution of water-waves</td><td>20</td></tr><tr><td>3.1. Determination of γ or μ from the boundary condition</td><td>20</td></tr><tr><td>3.2. Displacement of the free surface</td><td>24</td></tr><tr><td>3.3. Bernoulli equation and dipole formulation</td><td>26</td></tr><tr><td>3.4. Euler equation and vortex formulation</td><td>29</td></tr><tr><td>3.5. The deep-water case</td><td>32</td></tr><tr><td>4. Numerical discretization</td><td>33</td></tr><tr><td>4.1. Time integration</td><td>33</td></tr><tr><td>4.2. Shifted grids in space</td><td>34</td></tr><tr><td>4.3. Discretization and inversion of singular operator by Neumann series</td><td>34</td></tr><tr><td>5. Numerical results and convergence</td><td>37</td></tr><tr><td>5.1. Case 1: linear water waves</td><td>38</td></tr><tr><td>5.2. Case 2: solitary waves</td><td>39</td></tr><tr><td>5.3. Case 3: wave breaking</td><td>41</td></tr><tr><td>6. Comparison with regularization strategies</td><td>43</td></tr><tr><td>6.1. Fourier filtering for the vortex method</td><td>43</td></tr></table>

Date: December 12, 2023.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. The domain being considered consists of $\mathcal{D} = \mathbb{T}_L \times \mathbb{R}$, which is decomposed in $\mathcal{D} = \mathcal{D}_F \cup \mathcal{D}_A \cup \mathcal{D}_B \cup \Gamma_S \cup \Gamma_B$, where $\Gamma_S = \overline{\mathcal{D}_A} \cap \overline{\mathcal{D}_F}$, and $\Gamma_B = \overline{\mathcal{D}_F} \cap \overline{\mathcal{D}_B}$.</figcaption></figure>

<table><tr><td>6.2. Curve-offset method</td><td>44</td></tr><tr><td>7. Discussion</td><td>47</td></tr><tr><td>Acknowledgements</td><td>48</td></tr><tr><td>Appendix A. Cotangent kernel</td><td>48</td></tr><tr><td>Appendix B. Discrete operators for the dipole formulation</td><td>49</td></tr><tr><td>Appendix C. Discrete operators for the vortex formulation</td><td>51</td></tr><tr><td>Appendix D. Notations</td><td>54</td></tr><tr><td>References</td><td>55</td></tr></table>

# 1. INTRODUCTION

The study of water waves has a long mathematical history (Airy, Boussinesq, Cauchy, Kelvin, Laplace, Navier, Rayleigh, Saint-Venant, Stokes, to cite only a few). It has been studied in a variety of situations, probably the most complex of these being the wave breaking problem. What happens when a waves overturns raises significant mathematical difficulties. The water-air interface cannot be described as a graph any longer. A parametric description of the interface and the tracking of its Lagrangian evolution are needed.

We want to derive here a stable numerical strategy to solve for one-dimensional water waves (i.e. in a 2D domain, or a 3D domain assuming independence in one horizontal coordinate of space). We numerically approximate the free-surface Euler equations without introducing artificial regularizing parameters. This is particularly important in the case of loss of regularity of the interface, in order to study the possible formation of singularity (e.g. [7]).

1.1. **Problem formulation.** We consider a simple periodic domain $\mathcal{D} = \mathbb{T}_L \times \mathbb{R}$ see Fig. 1, and introduce two boundaries, $\Gamma_S$ the free surface water-air, and $\Gamma_B$ the bottom. The domain is thus decomposed in three subdomains, $\mathcal{D}_F$ the fluid domain, $\mathcal{D}_A$ the air domain, $\mathcal{D}_B$ below the bottom.
Since we want our mathematical approximation to be able to describe an interface which is not a graph (i.e. overturning of water in the context of a breaking wave) we need to be able to track it as a parametrized curve. Indeed, a description in the form $y = h(x)$ would develop a shock (discontinuity) as soon as the water starts to overturn.

Two cases will be considered, the single fluid problem, in which the air density is neglected, and the bi-fluid problem. In the latter case, the Euler equation needs to be considered both in the water ($D_F$) and in the air ($D_A$). At the water-bottom interface ($\Gamma_B$) the normal component of velocity needs to vanish (impermeability condition). Whereas at the water-air interface ($\Gamma_S$) two quantities need to be continuous across the interface: the normal velocity and the pressure. The latter, for an inviscid fluid, being equivalent to the continuity of the normal component of the stress tensor.

The velocity tangential to the interface is notably not continuous across $\Gamma_S$. This results in a localized distribution of vorticity along $\Gamma_S$ in the form of a vortex sheet.

The free surface $\Gamma_S$ is initially parametrized by arclength from left to right: $e \in [0, L_S] \mapsto z_S(e) = z_{S,1}(e) + iz_{S,2}(e)$. In the same way, the bottom $\Gamma_B$ is parametrized by $e \in [0, L_B] \mapsto z_B(e) = z_{B,1}(e) + iz_{B,2}(e)$.

It is useful to introduce the tangent vector $\tau_S = \tau_{S,1} + i\tau_{S,2} = |z_{S,e}(e)|^{-1}z_{S,e}$, which is pointing to the right and the normal $n_S = n_{S,1} + in_{S,2} = -\tau_{S,2} + i\tau_{S,1} = i\tau_S$ is pointing out of the fluid domain (where $z_{S,e}$ is the derivative in $e$). The same is done at the bottom with the normal now pointing in the fluid domain.

It should be stressed that the arclength is not preserved as the fluid surface $\Gamma_S$ evolves, it is thus important to consider $|z_e|(t, e)$.

We introduce on any vector field $\mathbf{u} = (u_1, u_2)$ the following three operations $\hat{\mathbf{u}} = u_1 - iu_2$, $(u_1, u_2)^\perp = (-u_2, u_1)$, and the curl operator $\mathrm{curl} \, \mathbf{u} = \partial_1 u_2 - \partial_2 u_1$. Finally, we also introduce for any vector $\mathbf{x} = (x_1, x_2)$ the complex notation $\mathbf{x} = x_1 + ix_2$.

## 1.2. Numerical strategies.
Interface evolution methods aim at capturing the time evolution of the water-air interface ($\Gamma_S$) using solely the knowledge of this vorticity distribution. This results in having to cope with singular integrals along the vortex sheet, but simplifies the problem numerically in that neither the water domain ($D_F$), nor the air domain ($D_A$), need to be meshed, as would be the case for example using a Finite-Element Method. Such approaches are closer in spirit to a Boundary-Element Method used in three-dimensional problems (e.g. [20, 21, 30, 31]).

A classical formulation of water waves is the celebrated Zakharov-Craig-Sulem description [37, 15]. While this description has proven extremely useful from a theoretical point of view (e.g. [2, 26, 18]), it raises difficulties from a numerical point of view and its adaptation to the framework of singular integral formulations could be the subject of a future work (see Remark 2.5 for a discussion).

We will consider instead two different approaches to compute the free-surface evolution. The first one is based on a so-called ‘vortex method’, the discontinuity of the tangential velocity at the air-water interface is modelled by a vortex sheet (with vorticity distribution $\gamma$). In this approach, a finite but large number of localized vortices is used to approximate a continuous vortex sheet. Such an approach has been shown to be efficient for Euler flows in the full space [19] and for exterior domains [4]. The interface is transported by
the resulting normal flow (as expressed in (3.3)). In the case of free surface flows, such as water waves, the vorticity distribution along the $Γ_S$ surface is however not preserved following a Lagrangian trajectory, and its evolution in time needs to be traced in the system describing the evolution of two Euler flows with two continuity conditions (see Equation (3.13)). This approach has been pioneered in [6] and further developed recently in [3]. It is also a useful description for mathematical proofs (e.g. [36, 14]).

The second approach is referred to as the ‘dipole layer’, where the velocity is now described by the jump in potential between the two fluids (measured by the dipole distribution $\mu$, related to $\gamma$ via $\gamma = \partial_e\mu$). The jump in tangential velocity then stems from expressing the velocity as a potential. The dipole layer evolution follows from the Bernoulli equation (and takes the form given in (3.10)), whereas the interface is again transported by the normal flow (as expressed in (3.4)). This approach is known as the ‘dipole method’ (it is equivalent to double layer potentials in potential theory). This approach was first investigated in [6] and then in [5, 7].

The former method is lighter to derive and offers the possibility to account for Euler flows with vorticity, whereas the latter involves more analytical work and assumes an irrotational flow. We will see however that the latter has better convergence properties for strongly non-linear configurations. In both cases, the spatial discretization in terms of singular integrals is known to converge toward the Euler equation [4], see also [8] for a study of the vortex method in the deep-water case.

The deep-water case is a trivial limit of the above description, in which the bottom vortex sheet (distributed on $Γ_B$) is sent to infinity. It is formally equivalent to simply suppressing it, or setting its vorticity to zero.

**1.3. Objectives of the present study.** Our aim is to construct a numerical scheme based on solid mathematical developments and free of smoothing or regularizations and which can later be used to guide theoretical understanding and further mathematical constructions on these problems. The goal of this article is thus to derive a formulation of these methods in the most general case (bi-fluid or single fluid, including possibly a non-flat bottom, including vorticity and mean currents).

The discrete approximation will be based on a relevant reformulation of the continuous problem. The first step is to introduce some quantities defined on the 1D interface, which then allows to recover the full fluid motion.

**Theorem 1.1.** Given $Γ_S, Γ_B ∈ C^{1,1}$, the following holds true.

For any $\mathbf{v} \in L^2(\Omega_F)$ such that $\text{div} \mathbf{v} = 0$, $\text{curl} \mathbf{v} \in L^\infty(\Omega_F)$, $\mathbf{v} \cdot n|_{Γ_B} = 0$, there exists $(\gamma_S, \gamma_B) \in C^{0,\alpha}$ and $(\mu_S, \mu_B) \in C^{1,\alpha}$ (for any $\alpha \in (0,1)$) such that $\mathbf{v}$ can be recovered uniquely in terms of $(\gamma_S, \gamma_B, \int_{Γ_B} \mathbf{v} \cdot \tau, \text{curl} v)$ or in terms of $(\mu_S, \mu_B, \int_{Γ_B} \mathbf{v} \cdot \tau, \text{curl} v)$.

This theorem will offer a base to turn a 2D problem into a 1D problem. It involves singular integral representation. This allows us to describe the free surface problem knowing only $(\gamma_S, \gamma_B)$ in the case of the vortex formulation, or $(\mu_S, \mu_B)$ in the case of the dipole formulation.

**Remark 1.2.** It will be showed in a second step that $\gamma_B$ (resp. $\mu_B$) can be determined uniquely in terms of $\gamma_S$ (resp. $\gamma_B$) and $(\int_{Γ_B} \mathbf{v} \cdot \tau, \text{curl} v)$.
The next step will then be to reformulate the water-waves equations into an evolution equation for $\gamma_S$ or $\mu_S$ and $z_S$.

**Theorem 1.3.** Let $z_S$ and $v$ be a regular solution of the water-waves equations, then the following holds:

(1) The single fluid equation without vorticity but with mean current and the bi-fluid model without vorticity and without mean current can be reformulated as explicit equations on $(\partial_t z_S, \partial_t \mu_S)$ which involves only $(z_S, \mu_S, \int_{\Gamma_B} v_0 \cdot \tau)$: see (3.4) and (3.10).

(2) The single fluid and bi-fluid equations with vorticity and with mean current can be reformulated as explicit equations on $(\partial_t z_S, \partial_t \gamma_S)$ which involves only $(z_S, \gamma_S, \int_{\Gamma_B} v_0 \cdot \tau, \text{curl } v)$: see (3.3) and (3.13).

The first item (1) corresponds to the dipole formulation, whereas the second (2) corresponds to the vortex formulation. Of course, to be complete, we need to add the transport equation for $\text{curl } v$ by the velocity given in terms of $(\gamma_S, \int_{\Gamma_B} v_0 \cdot \tau, \text{curl } v)$. It could appear strange to consider the dipole formulation which is less general than the vortex formulation, however we will see that it gives more stable numerical schemes. Surface tension is included in the derivation, but can also be omitted.

Both theorems above do not introduce any regularization and hold true for the continuous problem. The resulting expressions are thus non-trivial and, in a first reading, we advise to drop all terms associated to the density of air (i.e. associated to the bi-fluid formulation), uniform or localized vorticity and circulation (mean currents).

We then want to ensure that our numerical scheme converges in a realistic manner (i.e. for realistic parameters that can be achieved in practice) to the solution of the continuous problem. We introduce in this work a regularization-free approach to solve for the water-wave problem (i.e. without explicit filtering or any other regularization introducing extra parameters to the problem). To discretise singular integrals, we follow an approach introduced in [4] for which rigorous convergence proofs were provided in the case of a smooth boundary ($C^\infty$).

We verify conserved quantities at the discrete level. We illustrate on simple test cases the numerical convergence to the approximate solutions (e.g. Stokes waves, Green-Naghdi solitary waves). We also demonstrate stability and convergence of our numerical solution for the wave-breaking problem.

Finally, we investigate the effects of regularization strategies on the solution and illustrate numerically how they can yield irrelevant solutions.

**1.4. Plan of the paper.** In the next section, we introduce the singular integral representation, thus proving Theorem 1.1 by solving the corresponding elliptic equations.

In Section 3, we first prove Remark 1.2 and then reformulate the water-waves equations using this formalism, hence proving Theorem 1.3.

In these two sections, we have chosen to present a general derivation, using a minimal amount of simplifying assumptions. Simplifying assumptions are thus introduced as the derivation proceeds, only when they become necessary. We believe this highlights why
and where each hypothesis is needed. This justifies in particular the additional restrictions associated with the dipole formulation.

Section 4 presents our discretization strategy. Numerical results as well as convergence tests are presented in Section 5. The comparison with previously used regularization strategies (filtering and offsetting) is performed in Section 6. Finally in Section 7 we discuss potential applications and further development.

Finally, the Plemelj formulae, discrete expressions both for the vortex and the dipole method, as well as a list of notations are presented in appendices.

## 2. SINGULAR INTEGRALS REPRESENTATION AT FIXED TIME

In this section we will establish Theorem 1.1 by studying the elliptic problem related to inviscid flows.

### 2.1. Stream and potential functions.
The aim of this subsection is to express the velocity in both fluids in terms of stream-function and velocity potential. This involves the resolution of a div-curl problem in terms of the vorticity and the circulation.

#### 2.1.1. Resolution of the Laplace problem in $T_L \times \mathbb{R}$.
Even if we assume that the fluid is curl-free in $D_F$, the non-trivial boundary condition on $\Gamma_S$ will be interpreted as a vortex sheet in $D = T_L \times \mathbb{R}$. For this reason, we introduce the Green kernel in $D$

$$G(\mathbf{x}) = \frac{1}{4\pi} \ln \left( \cosh \frac{2\pi x_2}{L} - \cos \frac{2\pi x_1}{L} \right), \quad G(\mathbf{x}, \mathbf{y}) = G(\mathbf{x} - \mathbf{y}), \qquad (2.1)$$

and we recall the following result proved in [10]:

**Proposition 2.1.** For any $f \in L_c^\infty(D)$, every solution $\Psi$ of the following elliptic problem

$$\Delta\Psi = f, \quad \lim_{x_2 \to +\infty} \partial_2\Psi = - \lim_{x_2 \to -\infty} \partial_2\Psi, \quad |\Psi| \le C_1(|x_2| + 1)$$

can be written as

$$\Psi(\mathbf{x}) = \Psi[f](\mathbf{x}) = \int_D G(\mathbf{x}, \mathbf{y}) f(\mathbf{y}) \, d\mathbf{y} + C_2 \quad \text{where } C_1 \text{ and } C_2 \text{ are constants.} \qquad (2.2)$$

The relation between the Green kernel (2.1) in $T_L \times \mathbb{R}$ and the usual kernel $\frac{1}{2\pi} \ln |\mathbf{x} - \mathbf{y}|$ in $\mathbb{R}^2$ is formally derived in Appendix A.

From the explicit formula, it is easy to observe using Taylor expansions that

$$\begin{align}
\int_D G(\mathbf{x}, \mathbf{y}) f(\mathbf{y}) \, d\mathbf{y} &= \nonumber \\
&\quad \left(\frac{x_2}{2L} - \frac{\ln 2}{4\pi}\right) \int_D f(\mathbf{y}) \, d\mathbf{y} - \frac{1}{2L} \int_D y_2 f(\mathbf{y}) \, d\mathbf{y} + \mathcal{O}(e^{-|x_2|}) \quad \text{as } x_2 \to +\infty, \nonumber \\
\int_D G(\mathbf{x}, \mathbf{y}) f(\mathbf{y}) \, d\mathbf{y} &= \nonumber \\
&\quad \left(-\frac{x_2}{2L} - \frac{\ln 2}{4\pi}\right) \int_D f(\mathbf{y}) \, d\mathbf{y} + \frac{1}{2L} \int_D y_2 f(\mathbf{y}) \, d\mathbf{y} + \mathcal{O}(e^{-|x_2|}) \quad \text{as } x_2 \to -\infty. \tag{2.3}
\end{align}$$
Thus $\lim_{x_2 \to +\infty} \partial_2 \Psi = - \lim_{x_2 \to -\infty} \partial_2 \Psi$ is a necessary and sufficient condition to use the representation formula (2.2).

2.1.2. *Stream function and potential construction in the fluid domain.* We now apply this elliptic formalism to define the stream function for the vortex method. Let us consider the following elliptic problem on $\mathcal{D}_F$: for any functions $g$ with zero mean and $\omega$, we want to analyze a vector field $\mathbf{u}$ such that

$$ \mathrm{div} \, \mathbf{u} = 0 \text{ in } \mathcal{D}_F, \quad \mathrm{curl} \, \mathbf{u} = \omega \text{ in } \mathcal{D}_F, \quad \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B, \quad \mathbf{u} \cdot \mathbf{n} = g \text{ on } \Gamma_S, \quad (2.4) $$

which we want to extend in $\mathcal{D}$, in order to be able to use the above proposition.

In (2.4), the divergence free assumption stems from the incompressibility property and the third condition corresponds to the impermeability of the boundary at the bottom. By the Stokes formula, the fact $g$ has a zero mean is a necessary condition coming from these two assumptions.

Unfortunately, (2.4) has infinitely many solutions because of the harmonic vector field, also called the constant background current in [28], $\mathbf{H}$:

$$ \mathrm{div} \, \mathbf{H} = \mathrm{curl} \, \mathbf{H} = 0 \text{ in } \mathcal{D}_F, \quad \mathbf{H} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B \cup \Gamma_S $$

for instance when $\Gamma_B = \mathbb{T}_L \times \{-1\}$ and $\Gamma_S = \mathbb{T}_L \times \{0\}$, we have $\mathbf{H} = \mathbf{e}_1 \mathbf{1}_{\mathcal{D}_F}$.

In order to uniquely determine $\mathbf{u}$ from $\omega$, we need to prescribe the circulation either on the bottom or below the free surface, knowing that we have the following compatibility condition from the Stokes formula

$$ \int_{\Gamma_B} \mathbf{u} \cdot \tau \, d\sigma - \int_{\Gamma_S} \mathbf{u} \cdot \tau \, d\sigma = \int_{\mathcal{D}_F} \omega \, d\mathbf{x} \quad (2.5) $$

where the integrals are taken from left to right.

**Lemma 2.2.** For any $\omega \in L^\infty(\mathcal{D}_F)$, $g \in C^0(\Gamma_S)$ with zero mean value, and $\gamma \in \mathbb{R}$ given, there exists a unique $\mathbf{u} \in H^1(\mathcal{D}_F)$ such that

$$ \mathrm{div} \, \mathbf{u} = 0 \text{ in } \mathcal{D}_F, \quad \mathrm{curl} \, \mathbf{u} = \omega \text{ in } \mathcal{D}_F, \quad \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B, \\ \mathbf{u} \cdot \mathbf{n} = g \text{ on } \Gamma_S, \quad \int_{\Gamma_B} \mathbf{u} \cdot \tau \, d\sigma = \gamma. \quad (2.6) $$

Moreover, there exists a unique $\psi_F \in H^2(\mathcal{D}_F)$ up to an arbitrary constant, such that

$$ \mathbf{u} = \nabla^{\perp} \psi_F. $$

*Proof.* The proof of this lemma comes from standard ideas in elliptic theory and we give here only the main lines.

The existence comes from the existence of $\phi \in H^1(\mathcal{D}_F)$ solving the following Laplace problem with Neumann boundary condition (where the zero mean assumption is needed):

$\Delta\phi = 0$ in $\mathcal{D}_F$, $\partial_n\phi = -(\nabla^\perp\Psi + \alpha\mathbf{e}_1) \cdot \mathbf{n}$ on $\Gamma_B$, $\partial_n\phi = g - (\nabla^\perp\Psi + \alpha\mathbf{e}_1) \cdot \mathbf{n}$ on $\Gamma_S$ where $\Psi = \Psi[\omega]$ is defined in (2.2) and $\alpha \in \mathbb{R}$ such that $\int_{\Gamma_B} (\nabla^\perp\Psi + \alpha\mathbf{e}_1) \cdot \tau \, d\sigma = \gamma$. Indeed, it is then enough to set $\mathbf{u} = \nabla\phi + \nabla^\perp\Psi + \alpha\mathbf{e}_1$.
Having such a **u**, it is always possible to construct the stream function $\psi_F$, because $\text{div} \mathbf{u} = 0$ and $\int_{\Gamma_S} \mathbf{u} \cdot \mathbf{n} = \int_{\Gamma_B} \mathbf{u} \cdot \mathbf{n} = 0$ imply that $\int_{\Gamma} \mathbf{u}^{\perp} \cdot \mathbf{\tau} = 0$ for any closed loop $\Gamma$, which allows us to construct $\tilde{\psi}_F$, uniquely up to a constant.

The uniqueness can be deduced from the stream function. Indeed, denoting with tilde the difference of two solutions, we note that $\partial_{\tau} \tilde{\psi}_F = \tilde{\mathbf{u}} \cdot \mathbf{n} = 0$ implies that $\tilde{\psi}_F$ is constant on each component of the boundary and we conclude by integrating by parts:

$$ \int_{\mathcal{D}_F} |\tilde{\mathbf{u}}|^2 = \tilde{\psi}_F|_{\Gamma_B} \int_{\Gamma_B} \tilde{\mathbf{u}} \cdot \mathbf{\tau} \, ds - \tilde{\psi}_F|_{\Gamma_S} \int_{\Gamma_S} \tilde{\mathbf{u}} \cdot \mathbf{\tau} \, ds = 0. \quad \square $$

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

We should note that the conservation laws for the 2D Euler equations (including the circulation and the total vorticity) imply that $\int_{\Gamma_B} \mathbf{u} \cdot \mathbf{\tau}$ and $\int_{\Gamma_S} \mathbf{u} \cdot \mathbf{\tau}$ are both conserved quantities.

In the case of the dipole formulation, we need to write **u** as a gradient, which is possible only by subtracting the curl and the circulation parts. Of course, we could take advantage of the fact $\mathbf{u} - \nabla^{\perp} \Psi[\omega] - \frac{\gamma}{L} e_1$ is curl free with zero circulation, and can thus be written as a gradient in $\mathcal{D}_F$. Nevertheless this approach introduces additional difficulties.

For example, if we take into account the density of air, it will be crucial to properly define the air velocity field. However, for the single-fluid water-waves equations (in which the density of air is neglected), we are left with several possible choices, in particular stationary vector-fields could be used, thus simplifying the computation below. To underline where the properties of the vector fields are important, we stay general for now and we will introduce constraints as they become necessary.

Let us consider any $\mathbf{u}_{\omega, \gamma}$ such that

$$ \text{div} \mathbf{u}_{\omega, \gamma} = 0 \text{ in } \mathcal{D}_F, \quad \text{curl} \mathbf{u}_{\omega, \gamma} = \omega \text{ in } \mathcal{D}_F, \quad \int_{\Gamma_B} \mathbf{u}_{\omega, \gamma} \cdot \mathbf{n} \, ds = 0, \quad \int_{\Gamma_B} \mathbf{u}_{\omega, \gamma} \cdot \mathbf{\tau} \, ds = \gamma. \quad (2.7) $$

Finding **u** the solution of (2.6) is equivalent to look for $\mathbf{u}_R := \mathbf{u} - \mathbf{u}_{\omega, \gamma}$ which is $\text{div}$ and $\text{curl}$ free, without circulation and flux. The existence and uniqueness of $\mathbf{u}_R$ comes from Lemma 2.2. It can be written as the perpendicular gradient of a stream function, but also as the gradient of a potential function:

$$ \mathbf{u}_R = \nabla \phi_F = \nabla^{\perp} \tilde{\psi}_F \quad (2.8) $$

where $\phi_F$ and $\tilde{\psi}_F$ are uniquely determined, up to a constant. Here we use that $\mathbf{u}_R$ is curl free with zero circulation to state that $\int_{\Gamma} \mathbf{u} \cdot \mathbf{\tau} = 0$ for any closed loop $\Gamma$. Even if it may not seem natural to study $\tilde{\psi}_F$ instead of $\psi_F$, we will see below that $\tilde{\psi}_F$ is an interesting quantity to consider for the dipole formulation.

2.1.3. *Extension to the full domain.* Now that we have established Lemma 2.2 and (2.8), and in order to apply Proposition 2.1 to obtain a representation formula, we first need to extend continuously the potential $\phi_F$ or the stream functions $\psi_F$ or $\tilde{\psi}_F$. Extending the potential is related to the fluid charge method developed in [4]. This method is unfortunately not relevant for a free surface problem, see Remark 3.3.
We, therefore, prefer to extend the stream functions continuously. This is equivalent to assuming the continuity of the normal part of the velocity across the boundary. Such an extended vector field is divergence free in the whole domain $\mathcal{D}$, hence can be written using a stream function, and the boundaries can be interpreted as vortex sheets, corresponding to the jump in the tangential velocity.

At the bottom, we extend $\mathbf{u}$ in the simplest possible way, i.e. such that

$$ \operatorname{div} \mathbf{u} = \operatorname{curl} \mathbf{u} = 0 \text{ in } \mathcal{D}_B, \quad \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B, \quad \int_{\Gamma_B} \mathbf{u} \cdot \tau \, ds = 0, $$

which implies $\mathbf{u}|_{\mathcal{D}_B} = 0$.

This is equivalent to extending $\psi$ by the constant $\psi_F|_{\partial\Gamma_B}$ (and indeed, $\mathbf{u} \cdot \mathbf{n} = 0$ implies that $\psi_F$ is constant on $\Gamma_B$). In order to use Proposition 2.1, we have to extend the stream function in the air such that $u_2 \to 0$ as $x_2 \to +\infty$. Hence we extend it in the air $\mathcal{D}_A$ with the unique<fnref n="1" /> solution of

$$ \operatorname{div} \mathbf{u} = \operatorname{curl} \mathbf{u} = 0 \text{ in } \mathcal{D}_A, \quad \mathbf{u} \cdot \mathbf{n} = g \text{ on } \Gamma_A, $$

$$ |\mathbf{u}| \to 0 \text{ when } x_2 \to \infty, \quad \int_{\Gamma_A} \mathbf{u} \cdot \tau \, ds = 0. $$

Note that this equation is the physically relevant formulation if we are interested in the bi-fluid water-wave model, for which the continuity of normal velocity simply reflects that the two fluids are not mixing.

Note also that it could, in principle, be possible to add some vortices in the air. We should stress however that the circulation has to vanish at infinity in order to use Proposition 2.1, if not, we would have to change the extension below the bottom.

This extended vector field can be expressed as

$$ \mathbf{u} = \nabla^{\perp} \psi, \quad (2.9) $$

where $\psi$ is continuous in $\mathcal{D}$ and determined up to an arbitrary constant. This extension will be sufficient for the vortex formulation.

Regarding the derivation of the dipole formulation, we now have to extend $\mathbf{u}_R$. It will be convenient for the bottom condition to extend $\mathbf{u}_R$ by zero in $\mathcal{D}_B$ (see further down Remark 3.1). In order to achieve this, we must add the following assumption on $\mathbf{u}_{\omega,\gamma}$:

$$ \mathbf{u}_{\omega,\gamma} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B. \quad (2.10) $$

This assumption allows us to extend $\mathbf{u}_R$ in $\mathcal{D}_B$ by zero and $\tilde{\psi}$ by the constant $\tilde{\psi}_F|_{\partial\Gamma_B}$.

Again, for a compatibility at infinity, and in order to write $\mathbf{u}_R$ as a potential, we must extend $\mathbf{u}_R$ in the air by a vector field satisfying

$$ \operatorname{div} \mathbf{u}_R = \operatorname{curl} \mathbf{u}_R = 0 \text{ in } \mathcal{D}_A, \quad \mathbf{u}_R \cdot \mathbf{n} = g - \mathbf{u}_{\omega,\gamma} \cdot \mathbf{n} \text{ on } \Gamma_A, $$

$$ |\mathbf{u}_R| \to 0 \text{ when } x_2 \to \infty, \quad \int_{\Gamma_S} \mathbf{u}_R \cdot \tau \, ds = 0. $$

<footnote n="1">¹The existence and uniqueness can be proved via the double layer potential together with the Green kernel in $\mathbb{T}_L \times \mathbb{R}$ (see Proposition 2.1).</footnote>
From the above equations, we can write $\mathbf{u}_R = \nabla^{\perp} \tilde{\psi} = \nabla \phi$, where $\tilde{\psi}$ is continuous and uniquely defined up to an arbitrary constant. The potential $\phi$ jumps across $\Gamma_S$ and $\Gamma_B$, and we have complete freedom to choose independently the constants in each of the connected components: $D_B$, $D_F$ and $D_A$. These four constants will be determined below in order to be able to write $\psi$, $\tilde{\psi}$ and $\phi$ in the form of a singular integral by applying Proposition 2.1.

Using the uniqueness of the solution to the elliptic problem (2.6), it follows that

$$ \mathbf{u} = \mathbf{u}_{\omega,\gamma} + \nabla^{\perp} \tilde{\psi} = \mathbf{u}_{\omega,\gamma} + \nabla \phi \text{ in } D_F. \qquad (2.11) $$

Even if we have already properly defined the extension of $\tilde{\psi}$ in order to be able to use a Biot-Savart representation formula (2.2), we still need to discuss the expression of $\mathbf{u}_{\omega,\gamma}$ in $D_A$ to infer the value of $\mathbf{u}$ in the air. Such a discussion is postponed to the end of the next section.

**Remark 2.3.** We can apply the whole analysis of this paper to treat cases involving several submerged solids $S_k \in D_S$, simply constructing the harmonic vector such that

$$ \operatorname{div} \mathbf{H} = \operatorname{curl} \mathbf{H} = 0 \text{ in } D_F, \quad \mathbf{H} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B \cup_k \partial S_k, \quad \int_{\Gamma_B} \mathbf{H} \cdot \tau = \gamma_0, \quad \int_{\partial S_k} \mathbf{H} \cdot \tau = \gamma_k $$

where $\gamma_k$ is initially given. In the same way, if we are only interested by the single-fluid water-waves equations, we can simply construct $\mathbf{H}$ initially in $D_F \cup \Gamma_S \cup D_A$ and this problem can be solved in the dipole formulation. If $\omega = \gamma_0 = \gamma_1 = \gamma_k = 0$ for all $k$, the dipole formulation is possible for both the single-fluid and bi-fluid water-waves equations. Otherwise, we will need to use the vortex formulation where the inclusion of such solids is a minor modification of the numerical code. In the vortex formulation, we can even include the case where the solids are moving with a prescribed velocity and rotation by setting $\mathbf{H} \cdot \mathbf{n} = (\ell_i + r_i \mathbf{x}^{\perp}) \cdot \mathbf{n}$.

The case of immersed solids moving under the influence of the flow involves the computation of pressure forces at the boundary of the solid (see e.g. equation (5.4) in [3]). The case of a floating (partially immersed) solid would be even more challenging (see the recent developments in [9, 11, 22]).

**2.2. Potential and dipole formulae.** We now want to use Proposition 2.1 to express $\hat{\mathbf{u}} = u_1 - i u_2$ (in the vortex formulation), $\phi$ and $\tilde{\psi}$ (in the dipole formulation) as singular integrals (i.e. equations (2.14), (2.20), (2.22) below). In the previous subsection, we have constructed continuous $\psi$ or $\tilde{\psi}$ on $D$, where the perpendicular gradient is continuous on $D_B \cup D_F \cup D_A$, and its normal part is continuous across the interfaces $\Gamma_B$ and $\Gamma_S$. Extending $\mathbf{u}$ or $\mathbf{u}_R$ in this way ensures that $\operatorname{div} \mathbf{u} = 0$ in $D$, whereas the jump of the tangential part can be seen as a vortex sheet, namely

$$ \operatorname{curl} \mathbf{u} = \Delta \psi = \omega + |z_{S,e}|^{-1} \gamma_S \delta_{\Gamma_S} + |z_{B,e}|^{-1} \gamma_B \delta_{\Gamma_B} \text{ in } D $$

where

$$ \begin{align} \gamma_S(e) &:= |z_{S,e}(e)| \left[ \lim_{\mathbf{z} \in D_F \to \mathbf{z}_S(e)} \mathbf{u} - \lim_{z \in D_A \to z_S(e)} \mathbf{u} \right] \cdot \tau(e) \tag{2.12} \\ \gamma_B(e) &:= -|z_{B,e}(e)| \left[ \lim_{\mathbf{z} \in D_F \to \mathbf{z}_B(e)} \mathbf{u} \right] \cdot \tau(e) \nonumber \end{align} $$
are such that the mean value is

$$ \int (\omega + |z_{S,e}|^{-1}\gamma_S\delta_{\Gamma_S} + |z_{B,e}|^{-1}\gamma_B\delta_{\Gamma_B}) = 0, \quad (2.13) $$

see (2.5). These formulae and the following ones also hold replacing $\mathbf{u}$, $\omega$, $\gamma_S$, $\gamma_B$ by $\mathbf{u}_R$, $0$, $\tilde{\gamma}_S$, $\tilde{\gamma}_B$.

Proposition 2.1 implies that $\psi$ is determined up to a constant, which is fixed when we choose to represent<fnref n="2" /> it as follows

$$ \begin{aligned} \psi(\mathbf{x}) &= \int_{\Gamma_S} G(\mathbf{x}, \mathbf{y}) |z_{S,e}|^{-1} \gamma_S \, d\sigma(\mathbf{y}) + \int_{\Gamma_B} G(\mathbf{x}, \mathbf{y}) |z_{B,e}|^{-1} \gamma_B \, d\sigma(\mathbf{y}) + \int_{\mathcal{D}_F} G(\mathbf{x}, \mathbf{y}) \omega(\mathbf{y}) \, d\mathbf{y} \\ &= \int_0^{L_S} G(\mathbf{x}, \mathbf{z}_S(e)) \gamma_S(e) \, de + \int_0^{L_B} G(\mathbf{x}, \mathbf{z}_B(e)) \gamma_B(e) \, de + \int_{\mathcal{D}_F} G(\mathbf{x}, \mathbf{y}) \omega(\mathbf{y}) \, d\mathbf{y}. \end{aligned} $$

By the explicit formula of the Green kernel, we deduce from the previous formula the Biot-Savart law which yields the velocity $\mathbf{u} = \nabla^{\perp} \psi$ for all $x$ in $\mathcal{D}_F \cup \mathcal{D}_A \cup \mathcal{D}_S$:

$$ \begin{aligned} \hat{\mathbf{u}}(x) &= \int_0^{L_S} \gamma_S(e) \widehat{\nabla^{\perp} G}(\mathbf{x} - \mathbf{z}_S(e)) \, de + \int_0^{L_B} \gamma_B(e) \widehat{\nabla^{\perp} G}(\mathbf{x} - \mathbf{z}_B(e)) \, de \\ &\quad + \int_{\mathcal{D}_F} \widehat{\nabla^{\perp} G}(\mathbf{x}, \mathbf{y}) \omega(\mathbf{y}) \, d\mathbf{y} \\ &= \int_0^{L_S} \gamma_S(e) \frac{1}{2Li} \cot \left( \frac{x - z_S(e)}{L/\pi} \right) \, de + \int_0^{L_B} \gamma_B(e) \frac{1}{2Li} \cot \left( \frac{x - z_B(e)}{L/\pi} \right) \, de, \quad (2.14) \\ &\quad + \int_{\mathcal{D}_F} \frac{1}{2Li} \cot \left( \frac{x - y}{L/\pi} \right) \omega(\mathbf{y}) \, d\mathbf{y} \end{aligned} $$

because

$$ \widehat{\nabla^{\perp}} G(\mathbf{x}) = \frac{-\sinh \frac{x_2}{L/(2\pi)} - i \sin \frac{x_1}{L/(2\pi)}}{2L \left( \cosh \frac{x_2}{L/(2\pi)} - \cos \frac{x_1}{L/(2\pi)} \right)} = \frac{1}{2Li} \cot \left( \frac{x_1 + i x_2}{L/\pi} \right), \quad (2.15) $$

where we have used that $-\sinh b - i \sin a = -i(\sin a - \sin(i b)) = -2i \sin \frac{a-ib}{2} \cos \frac{a+ib}{2}$ and $\cosh b - \cos a = \cos ib - \cos a = 2 \sin \frac{a-ib}{2} \sin \frac{a+ib}{2}$. This formula with cotangent kernel is singular when $x$ goes to the boundary $\Gamma_S \cup \Gamma_B$. This is natural because it encodes the jump of the tangential part of the velocity. The limit formula, the so called Plemelj formulae, will play a crucial role in the sequel and are recalled in Appendix A. Another key tool presented in this Appendix is the following desingularization rule

$$ \begin{aligned} &pv \int \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) f(e') \, de' \\ &= \int \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) \frac{f(e')z_e(e) - f(e)z_e(e')}{z_e(e)} \, de', \quad (2.16) \end{aligned} $$

because it transforms a principal value integral into a classical integral of a smooth function. This exact relation will be systematically used in order to handle regular terms,

<footnote n="2">²Even if $\delta_{\Gamma}$ is not a bounded function, it belongs to $H^{-1}(\mathcal{D})$ where the well-posedness of elliptic problem is usually proven, and the formula can be rigorously established for $C^1$ curve, see [29, 25, 17].</footnote>
which can be integrated with greater accuracy, resulting in improved stability. It is worth stressing that this desingularization does not alter the accuracy of the scheme, as opposed to regularization technics. We would like to stress again that this periodic Biot-Savart law is formally related to the usual Biot-Savart law in $\mathbb{R}^2$: see Appendix A.

Theorem 1.1 is stated for bounded vorticities. Because the resulting equations will be later discretised, we restrict our attention below to a vorticity $\omega$ composed of a constant part $\omega_0 1_{\mathcal{D}_F}$ and a part that we approximate by a sum of Dirac masses $\sum_{j=1}^{N_v} \gamma_{v,j} \delta_{z_{v,j}(t)}$, see [19].

The velocity generated by the Dirac masses is simply $\frac{1}{2Li} \sum_{j=1}^{N_v} \gamma_{v,j} \cot \left( \frac{x-z_{v,j}}{L/\pi} \right)$. The velocity associated to the constant part can be simplified thanks to an integration by parts

$$ 
\begin{align*}
\int_{\mathcal{D}_F} \nabla^{\perp} G(\mathbf{x} - \mathbf{y}) \omega_0 \, d\mathbf{y} &= \omega_0 \left( - \int_{\mathcal{D}_F} (\nabla G)(\mathbf{x} - \mathbf{y}) \cdot \mathbf{e}_2(\mathbf{y}) \, d\mathbf{y} \right) \\
&= \omega_0 \left( \int_{\mathcal{D}_F} \nabla_y (G(\mathbf{x} - \mathbf{y})) \cdot \mathbf{e}_2(\mathbf{y}) \, d\mathbf{y} \right) \\
&\quad - \int_{\mathcal{D}_F} \nabla_y (G(\mathbf{x} - \mathbf{y})) \cdot \mathbf{e}_1(\mathbf{y}) \, d\mathbf{y} \\
&= \omega_0 \left( \int_{\partial \mathcal{D}_F} G(\mathbf{x} - \mathbf{y}) \mathbf{e}_2(\mathbf{y}) \cdot \tilde{\mathbf{n}}_F(\mathbf{y}) \, d\sigma(\mathbf{y}) \right) \\
&\quad - \int_{\partial \mathcal{D}_F} G(\mathbf{x} - \mathbf{y}) \mathbf{e}_1(\mathbf{y}) \cdot \tilde{\mathbf{n}}_F(\mathbf{y}) \, d\sigma(\mathbf{y}) \\
&= - \frac{\omega_0}{4\pi} \int_{\partial \mathcal{D}_F} \ln \left( \cosh \frac{x_2 - y_2}{L/(2\pi)} - \cos \frac{x_1 - y_1}{L/(2\pi)} \right) \tilde{\mathbf{n}}_F^{\perp}(\mathbf{y}) \, d\sigma(\mathbf{y})
\end{align*}
$$

where $\tilde{\mathbf{n}}_F$ is the unit normal vector outward to $\mathcal{D}_F$. This implies that

$$ 
\begin{align}
& \int_{\mathcal{D}_F} \frac{1}{2Li} \cot\left(\frac{x-y}{L/\pi}\right) \omega_0 \, d\mathbf{y} \nonumber \\
&= \frac{\omega_0}{4\pi} \int_0^{L_S} \ln\left(\cosh \operatorname{Im} \frac{x-z_S(e)}{L/(2\pi)} - \cos \operatorname{Re} \frac{x-z_S(e)}{L/(2\pi)}\right) \overline{z_{S,e}(e)} \, de \nonumber \\
&\quad - \frac{\omega_0}{4\pi} \int_0^{L_B} \ln\left(\cosh \operatorname{Im} \frac{x-z_B(e)}{L/(2\pi)} - \cos \operatorname{Re} \frac{x-z_B(e)}{L/(2\pi)}\right) \overline{z_{B,e}(e)} \, de, \tag{2.17}
\end{align}
$$

which is well defined and continuous in $\mathcal{D}$.

Let us note that we can compute $\mathbf{u}_{\omega,\gamma}$ for $\omega = \omega_0 + \sum \gamma_{v,j} \delta_{z_{v,j}}$ in the same way.

Therefore, we have a complete formula (2.14) which gives $\mathbf{u} = \nabla^{\perp} \psi$ in terms of $\omega, \gamma_S$ and $\gamma_B$, which will be used for the vortex formulation.

For the dipole formulation, we have, exactly in the same way,

$$ 
\begin{align}
\widehat{\mathbf{u}}_R(x) = & \int_0^{L_S} \tilde{\gamma}_S(e) \frac{1}{2Li} \cot \left( \frac{x - z_S(e)}{L/\pi} \right) \, de \nonumber \\
& + \int_0^{L_B} \tilde{\gamma}_B(e) \frac{1}{2Li} \cot \left( \frac{x - z_B(e)}{L/\pi} \right) \, de, \tag{2.18}
\end{align}
$$

where

$$ \tilde{\gamma}_S(e) := |z_{S,e}(e)| \left[ \lim_{\mathbf{z} \in \mathcal{D}_F \to \mathbf{z}_S(e)} \mathbf{u}_R - \lim_{\mathbf{z} \in \mathcal{D}_A \to \mathbf{z}_S(e)} \mathbf{u}_R \right] \cdot \tau(e) $$
$$ \tilde{\gamma}_B(e) := -|z_{B,e}(e)| \left[ \lim_{\mathbf{z} \in \mathcal{D}_F \to \mathbf{z}_B(e)} \mathbf{u}_R \right] \cdot \tau(e). $$

In the previous subsection, we have defined $\mathbf{u}_R$ and the extension such that $\mathbf{u}_R = \nabla\phi$ in $\mathcal{D}_B \cup \mathcal{D}_F \cap \mathcal{D}_A$ where we have the choice to fix one constant by connected component. As the mean values of $\tilde{\gamma}_S$ and $\tilde{\gamma}_B$ are zero, we know from the behavior at infinity (2.3) that $\nabla\phi = \mathbf{u}_R = \nabla^{\perp}\tilde{\psi}$ goes to zero exponentially fast when $x_2 \to \infty$. In order to control the boundary term in the following computation, we thus set the constant in $\mathcal{D}_A$ such that $\phi$ goes to zero at infinity. In the same way, we set the constant in $\mathcal{D}_B$ so that $\phi_B \to 0$ when $x_2 \to -\infty$. As $\phi$ is not continuous across the interfaces and we need the value from both side, we denote the restriction of $\phi$ in $\mathcal{D}_F$ (resp. in $\mathcal{D}_A$ and in $\mathcal{D}_B$) by $\phi_F$ (resp. by $\phi_A$ and $\phi_B$). For any $\mathbf{x} \in \mathcal{D}_F$, we compute

$$ \begin{aligned} \phi(\mathbf{x}) &= \langle \phi_F, \Delta G(\cdot - \mathbf{x}) \rangle \\ &= - \int_{\mathcal{D}_F} \nabla \phi_F(\mathbf{y}) \cdot \nabla G(\mathbf{y} - \mathbf{x}) \, d\mathbf{y} + \int_{\Gamma_S} \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\ &\quad - \int_{\Gamma_B} \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\ &= \int_{\Gamma_S} \left( \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_F(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}) \\ &\quad - \int_{\Gamma_B} \left( \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_F(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}), \end{aligned} $$

where we keep in mind that $\mathbf{n} = \tau^{\perp}$ is pointing outward on $\Gamma_S$ whereas it is pointing inward on $\Gamma_B$. As $\mathbf{u}_R \cdot \mathbf{n}$ is continuous, we have

$$ \phi(\mathbf{x}) = \int_{\Gamma_S} \left( \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_A(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}) \\ - \int_{\Gamma_B} \left( \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_B(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}). $$

As $\Delta G(\cdot - \mathbf{x}) = 0$ in $\mathcal{D}_A \cup \mathcal{D}_B$ (for $\mathbf{x} \in \mathcal{D}_F$), we can integrate by parts in the air and in the bottom domains as we did above in $\mathcal{D}_F$ to state

$$ \begin{aligned} & \int_{\Gamma_S} \left( \phi_A(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_A(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}) = 0 \\ & \int_{\Gamma_B} \left( \phi_B(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_B(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \right) d\sigma(\mathbf{y}) = 0 \end{aligned} $$

where we have used the fact that $G(\mathbf{x}) = \mathcal{O}(x_2)$ and $\nabla G(\mathbf{x}) = \mathcal{O}(1)$ at infinity. This implies

$$ \phi(\mathbf{x}) = \int_{\Gamma_S} (\phi_F(\mathbf{y}) - \phi_A(\mathbf{y})) \partial_n G(\mathbf{y} - \mathbf{x}) d\sigma(\mathbf{y}) - \int_{\Gamma_B} (\phi_F(\mathbf{y}) - \phi_B(\mathbf{y})) \partial_n G(\mathbf{y} - \mathbf{x}) d\sigma(\mathbf{y}). $$
Doing a similar computation for $\mathbf{x} \in \mathcal{D}_A$:

$$
\begin{align*}
\phi(\mathbf{x}) &= \langle \phi_A, \Delta G(\cdot - \mathbf{x}) \rangle = - \int_{\mathcal{D}_A} \nabla \phi_A(\mathbf{y}) \cdot \nabla G(\mathbf{y} - \mathbf{x}) \, d\mathbf{y} - \int_{\Gamma_S} \phi_A(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\
&= \int_{\Gamma_S} \partial_n \phi_A(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) - \int_{\Gamma_S} \phi_A(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\
&= \int_{\Gamma_S} \partial_n \phi_F(\mathbf{y}) G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) - \int_{\Gamma_S} \phi_A(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\
&= \int_{\Gamma_S} \phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) - \int_{\Gamma_S} \phi_A(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\
&\quad - \int_{\Gamma_B} (\phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_F(\mathbf{y}) G(\mathbf{y} - \mathbf{x})) \, d\sigma(\mathbf{y}) \\
&= \int_{\Gamma_S} (\phi_F(\mathbf{y}) - \phi_A(\mathbf{y})) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) \\
&\quad - \int_{\Gamma_B} (\phi_F(\mathbf{y}) \partial_n G(\mathbf{y} - \mathbf{x}) - \partial_n \phi_B(\mathbf{y}) G(\mathbf{y} - \mathbf{x})) \, d\sigma(\mathbf{y}) \\
&= \int_{\Gamma_S} (\phi_F(\mathbf{y}) - \phi_A(\mathbf{y})) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}) - \int_{\Gamma_B} (\phi_F(\mathbf{y}) - \phi_B(\mathbf{y})) \partial_n G(\mathbf{y} - \mathbf{x}) \, d\sigma(\mathbf{y}),
\end{align*}
$$

we notice that this formula holds true in $\mathcal{D}_B \cup \mathcal{D}_F \cup \mathcal{D}_A$. So, we are computing now $\partial_n G(\mathbf{y} - \mathbf{x})$

$$
\begin{align*}
\nabla G(\mathbf{y} - \mathbf{x}) \cdot \mathbf{n}(\mathbf{y}) \, d\sigma(\mathbf{y}) &= \frac{1}{2L \left( \cosh \frac{z_2(e)-x_2}{L/(2\pi)} - \cos \frac{z_1(e)-x_1}{L/(2\pi)} \right)} \begin{pmatrix} \sin \frac{z_1(e)-x_1}{L/(2\pi)} \\ \sinh \frac{z_2(e)-x_2}{L/(2\pi)} \end{pmatrix} \cdot \begin{pmatrix} -z_{2,e}(e) \\ z_{1,e}(e) \end{pmatrix} \, de \\
&= - \operatorname{Re} \left[ \frac{-\sinh \frac{z_2(e)-x_2}{L/(2\pi)} - i \sinh \frac{z_1(e)-x_1}{L/(2\pi)}}{2L \left( \cosh \frac{z_2(e)-x_2}{L/(2\pi)} - \cos \frac{z_1(e)-x_1}{L/(2\pi)} \right)} z_e(e) \right] \, de \\
&= - \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z(e)-x}{L/\pi} \right) z_e(e) \right] \, de
\end{align*}
$$

so, setting

$$
\mu_S(e) = (\phi_F - \phi_A)(z_S(e)), \quad \mu_B(e) = (\phi_B - \phi_F)(z_B(e)), \qquad (2.19)
$$

we finally get for any $\mathbf{x} \in \mathcal{D}_F \cup \mathcal{D}_A \cup \mathcal{D}_B$

$$
\begin{align}
\phi(\mathbf{x}) &= - \int_0^{L_S} \mu_S(e) \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - x}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e \nonumber \\
&\quad - \int_0^{L_B} \mu_B(e) \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - x}{L/\pi} \right) z_{B,e}(e) \right] \mathrm{d}e \nonumber \\
&= \int_0^{L_S} \mu_S(e) \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_S(e)}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e \tag{2.20}
\end{align}
$$
$$ + \int_{0}^{L_B} \mu_B(e) \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_B(e)}{L/\pi} \right) z_{B,e}(e) \right] \mathrm{d}e. $$

We note here that we did not provide any restriction on the constant for $\phi_F$ so the previous formula holds true if we change $\phi_F$ (so $\mu_S$ and $\mu_B$) by a constant. It is therefore possible to fix initially this constant in such a way that

$$ \int_{0}^{L_S} \mu_{S,0}(e) \mathrm{d}e = 0. \qquad (2.21) $$

This condition is not conserved in time.

Let us also note that with our extension and Assumption (2.10), we have $\phi_B = 0$ in $D_B$.

It is also possible to derive the stream function $\tilde{\psi}$ from $\mu_S$ and $\mu_B$. To do this, we first remark that

$$ \mathbf{u}_R \cdot \tau = \nabla \phi \cdot \tau = |z_e|^{-1} \partial_e (\phi(z)), $$

hence

$$ \tilde{\gamma}_S(e) = \partial_e \left[ \phi_F(z_S(e)) - \phi_A(z_S(e)) \right] = \partial_e \mu_S(e) \quad \text{and} \quad \tilde{\gamma}_B(e) = -\partial_e \phi_F(z_B(e)) = \partial_e \mu_B(e), $$

and then for any constants $C_S, C_B \in \mathbb{R}$

$$ \begin{aligned} \tilde{\psi}(x) &= \int_{0}^{L_S} \partial_e(\mu_S(e) + C_S)G(\mathbf{z}_S(e) - \mathbf{x}) \mathrm{d}e + \int_{0}^{L_B} \partial_e(\mu_B(e) + C_B)G(\mathbf{z}_B(e) - \mathbf{x}) \mathrm{d}e \\ &= -\int_{0}^{L_S} (\mu_S(e) + C_S)\partial_e\left(G(\mathbf{z}_S(e) - \mathbf{x})\right) \mathrm{d}e - \int_{0}^{L_B} (\mu_B(e) + C_B)\partial_e\left(G(\mathbf{z}_B(e) - \mathbf{x})\right) \mathrm{d}e. \end{aligned} $$

So we need to compute

$$ \begin{aligned} \nabla G(\mathbf{z}(e) - \mathbf{x}) \cdot \begin{pmatrix} z_{e,1} \\ z_{e,2} \end{pmatrix} &= \frac{1}{2L \left( \cosh \frac{z_2(e)-x_2}{L/(2\pi)} - \cos \frac{z_1(e)-x_1}{L/(2\pi)} \right)} \begin{pmatrix} \sin \frac{z_1(e)-x_1}{L/(2\pi)} \\ \sinh \frac{z_2(e)-x_2}{L/(2\pi)} \end{pmatrix} \begin{pmatrix} z_{e,1}(e) \\ z_{e,2}(e) \end{pmatrix} \\ &= - \operatorname{Im} \left[ \frac{-\sinh \frac{z_2(e)-x_2}{L/(2\pi)} - i \sinh \frac{z_1(e)-x_1}{L/(2\pi)}}{2L \left( \cosh \frac{z_2(e)-x_2}{L/(2\pi)} - \cos \frac{z_1(e)-x_1}{L/(2\pi)} \right)} z_e(e) \right] \\ &= - \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{z(e)-x}{L/\pi} \right) z_e(e) \right] \end{aligned} $$

and we finally get

$$ \tilde{\psi}(x) = \int_{0}^{L_S} (\mu_S(e) + C_S) \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - x}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e \\ + \int_{0}^{L_B} (\mu_B(e) + C_B) \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - x}{L/\pi} \right) z_{B,e}(e) \right] \mathrm{d}e $$
$$ = - \int_{0}^{L_S} (\mu_S(e) + C_S) \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_S(e)}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e \quad (2.22) \\ - \int_{0}^{L_B} (\mu_B(e) + C_B) \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_B(e)}{L/\pi} \right) z_{B,e}(e) \right] \mathrm{d}e. $$

As this formula is valid for any values of $C_B$ and $C_S$, it holds true for $C_S = C_B = 0$ and besides

$$ \int_{0}^{L_S} \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_S(e)}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e = \int_{0}^{L_B} \operatorname{Im} \left[ \frac{1}{2Li} \cot \left( \frac{x - z_B(e)}{L/\pi} \right) z_{B,e}(e) \right] \mathrm{d}e = 0. $$

For Section 3.3, it will be convenient to introduce the quantity

$$ \Phi_S(e) = (\phi_F + \phi_A)(z_S(e)), $$

which is complementary to $\mu_S$, and which can be expressed thanks to the formula giving $\phi$ and the limit formula (see Appendix A)

$$ \Phi_S(e) = \int_0^{L_S} \mu_S(e') \operatorname{Re} \left[ \frac{1}{Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) z_{S,e}(e') \right] \mathrm{d}e' \\ + \int_0^{L_B} \mu_B(e') \operatorname{Re} \left[ \frac{1}{Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) z_{B,e}(e') \right] \mathrm{d}e'. $$

Computing the limit for $e' \to e$, we note that the first integral is a classical integral of a continuous function, where the extension for $e' = e$ is

$$ -\mu_S(e) \operatorname{Re} \left[ \frac{z_{S,ee}(e)}{2\pi i z_{S,e}(e)} \right]. $$

Even if this integral could be well approximated by Riemann sum for smooth fluid surface, it occurs that the following formula will be convenient to get non singular integrals, taking advantage of the desingularization (2.16)

$$ \Phi_S(e) = \int_0^{L_S} (\mu_S(e') - \mu_S(e)) \operatorname{Re} \left[ \frac{1}{Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) z_{S,e}(e') \right] \mathrm{d}e' \\ + \int_0^{L_B} \mu_B(e') \operatorname{Re} \left[ \frac{1}{Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) z_{B,e}(e') \right] \mathrm{d}e' \quad (2.23) $$

which is then extended for $e = e'$ by zero.

We conclude this section with one last compatibility condition, which is not used in this article but will be used in a forthcoming article.

**Remark 2.4.** The function $z \mapsto \phi - i\tilde{\psi}$ is harmonic in $D_B \cup D_F \cup D_A$, hence the integrals along two curves going from left to right are the same if both curves are included in the
same connected component. With the limit at infinity, it is clear that

$$ \int_0^{L_S} (\phi_A(z_S(e)) - i\tilde{\psi}(z_S(e))) z_{S,e}(e) \, de = -iL \lim_{x_2 \to +\infty} \tilde{\psi}. $$

whereas

$$ \int_0^{L_B} (\phi_B(z_B(e)) - i\tilde{\psi}(z_B(e))) z_{B,e}(e) \, de = -iL \lim_{x_2 \to -\infty} \tilde{\psi} = iL \lim_{x_2 \to +\infty} \tilde{\psi}. $$

Inside the fluid we have

$$ \int_0^{L_S} (\phi_F(z_S(e)) - i\tilde{\psi}(z_S(e))) z_{S,e}(e) \, de = \int_0^{L_B} (\phi_F(z_B(e)) - i\tilde{\psi}(z_B(e))) z_{B,e}(e) \, de. $$

By continuity of the stream function, we get

$$ \begin{aligned} \int_0^{L_S} \mu_S(e) z_{S,e}(e) \, de &= \int_0^{L_S} (\phi_F - \phi_A)(z_S(e)) z_{S,e}(e) \, de \\ &= \int_0^{L_S} (\phi_F(z_S(e)) - i\tilde{\psi}(z_S(e))) z_{S,e}(e) \, de \\ &\quad - \int_0^{L_S} (\phi_A(z_S(e)) - i\tilde{\psi}(z_S(e))) z_{S,e}(e) \, de \\ &= \int_0^{L_B} (\phi_F(z_B(e)) - i\tilde{\psi}(z_B(e))) z_{B,e}(e) \, de + iL \lim_{x_2 \to +\infty} \tilde{\psi} \\ &= \int_0^{L_B} (\phi_F - \phi_B)(z_B(e)) z_{B,e}(e) \, de + 2iL \lim_{x_2 \to +\infty} \tilde{\psi} \\ &= - \int_0^{L_B} \mu_B(e) z_{B,e}(e) \, de + 2iL \lim_{x_2 \to +\infty} \tilde{\psi} \end{aligned} $$

we thus have for all time

$$ \int_0^{L_S} \mu_S(e) \operatorname{Re}[z_{S,e}(e)] \, de = - \int_0^{L_B} \mu_B(e) \operatorname{Re}[z_{B,e}(e)] \, de. \quad (2.24) $$

*Remark 2.5.* The celebrated Dirichlet to Neumann operator in the Zakharov-Craig-Sulem formulation [37, 15] is very close to the dipole derivation. For $\varphi \in H^{1/2}(\Gamma_S)$ given, the principle is indeed to find $u_R = \nabla\phi_F = \nabla^\perp\tilde{\psi}$ such that

$$ \Delta\phi_F = 0 \text{ in } \mathcal{D}_F, \partial_n\phi_F = 0 \text{ on } \Gamma_B, \phi_F = \varphi \text{ on } \Gamma_S. $$

Extending as we did $\tilde{\psi}$ by continuity and defining $\phi$, we can represent $\phi$ through the singular representation formulation (2.20). Therefore, we should first find uniquely $\mu_S$ and $\mu_B$ such that $\phi_B = 0$ on $\Gamma_B$ and $\phi_F = \varphi$ on $\Gamma_S$ thanks to the limit formulae of Appendix A (see (3.2) for this kind of application). With $(\mu_S, \mu_B)$ found, we differentiate in order to get $(\tilde{\gamma}_S, \tilde{\gamma}_B)$ which allow us to construct $u_R$ (2.18), hence $\partial_n\phi_F|_{\Gamma_S}$ again with the limit formulae. This ends the definition of the Dirichlet to Neumann operator $\varphi \mapsto \partial_n\phi_F|_{\Gamma_S}$.
2.3. **Discussion about $\mathbf{u}_{\omega,\gamma}$ for the dipole formulation.** The simplest case that we will study in details is the case where $\omega = \gamma = 0$ where we choose of course $\mathbf{u}_{\omega,\gamma} = 0$. Then $\mathbf{u} = \nabla\phi_F = \nabla^{\perp}\tilde{\psi}_F$ is naturally defined in the full domain $\mathcal{D}$. In this easy case, we can consider both the bi-fluid water-wave equation, in which the air is assumed to be an incompressible fluid with a non-zero density, or the single-fluid water-waves equations, where we neglect the density of the air in $\mathcal{D}_A$.

When we have some vorticity or background current, we need to discuss the expression of $\mathbf{u}_{\omega,\gamma}$ in $\mathcal{D}_A$ to infer the value of $\mathbf{u}$ in the air. There are essentially two natural options:

* either to have an explicit formula for $\mathbf{u}_{\omega,\gamma}$, or at least assume it is independent of time;

* or to extend by zero.

The choice depends on whether we need the physical air velocity, i.e. if we consider the single fluid or the bi-fluid water-wave equation.

In the latter case (single-fluid), we do not need to know the velocity in the air, and we can simply set

$$ \mathbf{u}_{\omega,\gamma} = \frac{\gamma}{L} \mathbf{e}_1 \mathbf{1}_{\mathcal{D}_F} \quad \text{if the bottom is flat and } \omega = 0. $$

Then we cannot say that $\mathbf{u}_{\omega,\gamma} + \nabla^{\perp}\tilde{\psi}$ defines the velocity in the air, because the normal part of the velocity is not continuous. A natural idea would then be to set

$$ \mathbf{u}_{\omega,\gamma} = \frac{\gamma}{L} \mathbf{e}_1 \chi(x_2) \quad \text{if the bottom is flat and } \omega = 0. $$

If we choose $\chi(x_2) = 1$ for all $x_2$, this implies that a non-physical circulation is present in the air, which is equal to the circulation in the water. Alternatively if $\chi(x_2)$ is chosen to decay smoothly from 1 near the interface to 0 at infinity, this implies a strange, also non-physical, vorticity in the air curl $\mathbf{u} = -\frac{\gamma}{L}\chi'(x_2)$. Both cases do not correspond to the actual air velocity. Hence, in the limiting case of vanishing air density, we can use this simpler expression for the velocity $\mathbf{u}_{\omega,\gamma}$. The air velocity can, however, not be reconstructed in that case (as it does not influence the interface evolution). Therefore, we will consider later the case of the single-fluid water-waves equation without vorticity but with background current, constructing the time independent $\mathbf{H}$ solving

$$ \mathrm{div} \, \mathbf{H} = \mathrm{curl} \, \mathbf{H} = 0 \text{ in } \mathcal{D}_F \cup \Gamma_S \cup \mathcal{D}_A, \quad \mathbf{H} \cdot \mathbf{n} = 0 \text{ on } \Gamma_B, \quad \int_{\Gamma_B} \mathbf{H} \cdot \mathbf{\tau} = \gamma $$

and setting $\mathbf{u}_{\omega,\gamma} = \mathbf{H}$. Namely, we set

$$ \widehat{\mathbf{u}}_{\omega,\gamma}(x) = \widehat{\mathbf{H}}(x) = \int_{0}^{L_B} \gamma_{B,H}(e) \frac{1}{2Li} \cot\left(\frac{x - z_B(e)}{L/\pi}\right) \mathrm{d}e $$

where $\gamma_{B,H}$ is the unique<fnref n="3" /> solution of

$$ \mathrm{pv} \int_{0}^{L_B} \gamma_{B,H}(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) \right] \mathrm{d}e' = 0, $$

<footnote n="3">³See Section 3.1.2.</footnote>
with

$$ \int_{0}^{L_B} \gamma_{B,H}(e') \, de' = -2\gamma. $$

Indeed, $\widehat{\mathbf{u}}_{\omega,\gamma}$ constructed in this way has a circulation $\gamma$ in $D_F \cup \Gamma_S \cup D_A$ and $-\gamma$ in $D_B$, which is compatible with the limit behavior of the stream function associated to a vorticity which has a non vanishing mean value (see Proposition 2.1). For the flat bottom, we recover $\mathbf{H} = \frac{\gamma}{L}\mathbf{e}_1$ in $D_F \cup \Gamma_S \cup D_A$ and $\mathbf{H} = -\frac{\gamma}{L}\mathbf{e}_1$ in $D_S$.

In the case of the single-fluid water-waves equations with a flat bottom $\Gamma_B = \mathbb{T}_L \times \{-h_0\}$, we could think to the simple formula coming from the image method:

$$ \widehat{\mathbf{u}}_{\omega,\gamma}(x) = \frac{\gamma}{L} + \int_{D} \frac{1}{2Li} \cot\left(\frac{x-y}{L/\pi}\right) (\omega 1_{D_F} + \tilde{\omega} 1_{D_B})(\mathbf{y}) \, d\mathbf{y} $$

where $\tilde{\omega}(x_1, x_2) := -\omega(x_1, -x_2 - 2h_0)$. Hence, we have after two integrations by parts

$$ \begin{aligned} \widehat{\mathbf{u}}_{\omega,\gamma}(x) = & \frac{\gamma}{L} + \frac{1}{2Li} \sum_{j=1}^{N_v} \gamma_{v,j} \left( \cot\left(\frac{x-z_{v,j}}{L/\pi}\right) - \cot\left(\frac{x-\overline{z_{v,j}}+2ih_0}{L/\pi}\right) \right) \\ & + \frac{\omega_0}{4\pi} \int_0^{L_S} \ln\left(\cosh \operatorname{Im} \frac{x-z_S(e)}{L/(2\pi)} - \cos \operatorname{Re} \frac{x-z_S(e)}{L/(2\pi)}\right) \overline{z_{S,e}(e)} \, de \\ & - \frac{\omega_0}{2\pi} \int_0^{L_B} \ln\left(\cosh \operatorname{Im} \frac{x-z_B(e)}{L/(2\pi)} - \cos \operatorname{Re} \frac{x-z_B(e)}{L/(2\pi)}\right) \overline{z_{B,e}(e)} \, de \\ & + \frac{\omega_0}{4\pi} \int_0^{L_S} \ln\left(\cosh \operatorname{Im} \frac{x-\overline{z_S(e)}+2ih_0}{L/(2\pi)} - \cos \operatorname{Re} \frac{x-\overline{z_S(e)}+2ih_0}{L/(2\pi)}\right) z_{S,e}(e) \, de. \end{aligned} $$

Unfortunately, we will observe later (see Section 3.3) that this approach is unpractical.

In the first case (bi-fluid formulation), if we want to extend $\mathbf{u}$ in such a way that the normal component of the velocity is continuous and $\operatorname{div} \mathbf{u} = \operatorname{curl} \mathbf{u} = 0$ in $D_A$, we then have to solve at any time an elliptic problem in $D_A$ to extend the flow $\mathbf{u}_{\omega,\gamma}$ in the correct way. Alternatively, we could prefer to extend $\mathbf{u}_{\omega,\gamma}$ by zero. In this case, we would need to add the following condition

$$ \mathbf{u}_{\omega,\gamma} \cdot \mathbf{n} = 0 \text{ on } \Gamma_S. \qquad (2.25) $$

and solve at any time the elliptic problem (2.7) in $D_F$ with (2.10) and (2.25). Namely, we set

$$ \begin{aligned} \widehat{\mathbf{u}}_{\omega,\gamma}(x) = & \int_{0}^{L_S} \gamma_{S,\omega,\gamma}(e) \frac{1}{2Li} \cot\left(\frac{x-z_S(e)}{L/\pi}\right) \, de + \int_{0}^{L_B} \gamma_{B,\omega,\gamma}(e) \frac{1}{2Li} \cot\left(\frac{x-z_B(e)}{L/\pi}\right) \, de \\ & + \int_{D_F} \frac{1}{2Li} \cot\left(\frac{x-y}{L/\pi}\right) \omega(\mathbf{y}) \, d\mathbf{y} \end{aligned} $$

where $(\gamma_{S,\omega,\gamma}, \gamma_{B,\omega,\gamma})$ is properly constructed at each time. It will also be observed in Section 3.3 that this approach is too complicated. In the case of a bi-fluid formulation with circulation, or with internal vorticity, the vortex method will be preferred.

The dipole formulation can however be considered only in two cases: the case of a bi-fluid formulation in the absence of both internal vorticity and circulation or the case of a single fluid formulation in the absence of internal vorticity.
This concludes the proof of Theorem 1.1 since the solution of an elliptic problem in a $C^{1,1}$ domain belongs to $\cap_{p>2} W^{2,p}$ when $\omega \in L^\infty$. This justifies that the velocity and then $\gamma$ is $C^{0,\alpha}$ for all $\alpha \in (0,1)$. We recall that $\gamma_S = \partial\mu_S/\partial e$, which yields higher regularity on $\mu_S$.

Note that $\gamma_S \in C^{0,\alpha}$ is enough to be able to reformulate principal value integrals in the form of classical integrals (2.16).

### 3. EVOLUTION OF WATER-WAVES

The bottom $z_B$ and the constant part of the vorticity $\omega_0$ are initially given. At any time, for a given $(z_{v,j})_{j=1,...,N_v}$ and $z_S$, we have established in Section 2.1 the existence of $(\gamma_S, \gamma_B)$ or $(\mu_S, \mu_B)$ from $g$. Conversely, from $\gamma_S$ or $\mu_S$, we will first show that there is a unique $\gamma_B$ or $\mu_B$ satisfying the boundary conditions at the bottom, and hence proving Remark 1.2. We will thus use Section 2.2 to get the velocity everywhere, and then deduce the displacements of the point vortex and the free surface: $\partial_t z_{v,j}$ and $\partial_t z_S$. The last step is to use the Euler or the Bernoulli equations to determine $\partial_t \gamma_S$ or $\partial_t \mu_S$, namely proving Theorem 1.3.

Therefore, if we know $g$ initially, we can construct $(\gamma_{S,0}, \gamma_{B,0})$ or $(\tilde{\gamma}_{S,0}, \tilde{\gamma}_{B,0})$ such that the corresponding velocity (2.14) or (2.18) verifies the correct boundary conditions. From $\tilde{\gamma}_{S,0}$ we will construct $\mu_{S,0}$ as the primitive of $\tilde{\gamma}_{S,0}$ with zero mean. For $t > 0$, the main numerical strategy can be summarized as

* for the vortex method:

$$ (z_S, (z_{v,j})_j, \gamma_S) \mapsto (z_S, (z_{v,j})_j, \gamma_S, \gamma_B) \mapsto (\partial_t z_S, (\partial_t z_{v,j})_j, \partial_t \gamma_S); $$

* for the dipole method:

$$ (z_S, (z_{v,j})_j, \mu_S) \mapsto (z_S, (z_{v,j})_j, \mu_S, \mu_B) \mapsto (\partial_t z_S, (\partial_t z_{v,j})_j, \partial_t \mu_S). $$

### 3.1. Determination of $\gamma$ or $\mu$ from the boundary condition.
The quantities $z_B$, $\gamma$, $\omega_0$, $(\gamma_{v,j})_j$, $g$, are given by the initial conditions, and we want to solve

* $(z_{S,0}, (z_{v,j,0})_j, g_0) \mapsto \gamma_{S,0}$ for the initial setting in the vortex formulation;

* $(z_S, (z_{v,j})_j, \gamma_S) \mapsto \gamma_B$ for every time step in the vortex formulation;

* $(z_{S,0}, u_{\omega,\gamma,0}, g_0) \mapsto \tilde{\gamma}_{S,0} \mapsto \mu_{S,0}$ for the initial setting in the dipole formulation;

* $(z_S, u_{\omega,\gamma}, \mu_S) \mapsto \mu_B$ for every time step in the dipole formulation.

#### 3.1.1. Initial $\gamma_{S,0}$ for the vortex formulation.
In many situation, such as solitary waves, $z_B$, $z_S$, $g = u_F \cdot n|_{\Gamma_S}$, $\gamma$, $\omega_0$ and $(\gamma_{v,j}, z_{v,j})_{j=1,...,N_v}$ are known initially. By uniqueness of the elliptic problem (see (2.6) with our extension (2.9)), we know that there exists a unique pair $(\gamma_S, \gamma_B)$ such that the normal velocity $u \cdot n = - \operatorname{Im}(\hat{u}_{|z_S,e}^{z_{S,e}})$ verifies the proper boundary condition on the free surface. Using (2.14),

$$ \mathrm{pv} \int_{0}^{L_S} \gamma_S(e') \operatorname{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] de' $$
$$ + \int_{0}^{L_B} \gamma_B(e') \operatorname{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \right] de' = \operatorname{RHS}_{V0,S}(e), $$

where thanks to (2.17)

RHS$_{V0,S}$(e)

$$ \begin{aligned} &= -g|z_S(e)| - \int_{\mathcal{D}_F} \operatorname{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - y}{L/\pi} \right) \right] \omega_0 \, dy \\ &= -g|z_S(e)| - \sum_{j=1}^{N_v} \gamma_{v,j} \operatorname{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_{v,j}}{L/\pi} \right) \right] \\ &\quad - \frac{\omega_0}{4\pi} \int_0^{L_S} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_S(e')}{L/(2\pi)} \right) \operatorname{Im} \left[ z_{S,e}(e) \overline{z_{S,e}(e')} \right] de' \\ &\quad + \frac{\omega_0}{4\pi} \int_0^{L_B} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_B(e')}{L/(2\pi)} \right) \operatorname{Im} \left[ z_{S,e}(e) \overline{z_{B,e}(e')} \right] de'; \end{aligned} $$

on the bottom:

$$ \int_0^{L_S} \gamma_S(e') \operatorname{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\ + \operatorname{pv} \int_0^{L_B} \gamma_B(e') \operatorname{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) \right] de' = \operatorname{RHS}_{V0,B}(e) $$

where

RHS$_{V0,B}$(e)

$$ \begin{aligned} &= - \int_{\mathcal{D}_F} \operatorname{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - y}{L/\pi} \right) \right] \omega_0 \, dy \\ &= - \sum_{j=1}^{N_v} \gamma_{v,j} \operatorname{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_{v,j}}{L/\pi} \right) \right] \\ &\quad - \frac{\omega_0}{4\pi} \int_0^{L_S} \ln \left( \cosh \operatorname{Im} \frac{z_B(e) - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_B(e) - z_S(e')}{L/(2\pi)} \right) \operatorname{Im} \left[ z_{B,e}(e) \overline{z_{S,e}(e')} \right] de' \\ &\quad + \frac{\omega_0}{4\pi} \int_0^{L_B} \ln \left( \cosh \operatorname{Im} \frac{z_B(e) - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_B(e) - z_B(e')}{L/(2\pi)} \right) \operatorname{Im} \left[ z_{B,e}(e) \overline{z_{B,e}(e')} \right] de', \end{aligned} $$

together with the circulation assumptions:

$$ \int_{0}^{L_S} \gamma_S(e') de' = \gamma - \omega_0 |\mathcal{D}_F| - \sum_{j} \gamma_{v,j} $$
$$ \int_{0}^{L_B} \gamma_B(e') de' = -\gamma. $$

The existence and uniqueness of a solution is related to the operator B in [4], and Section 4.3 will detail how these integrals can be discretized, ensuring that the resulting matrices are invertible.

### 3.1.2. Time dependent $\gamma_B$ for the vortex formulation.
For any given time, knowing $z_B, \gamma$, $\omega_0$, $(\gamma_{v,j})_{j=1,...,N_v}$ from the initial conditions, we need to construct $\gamma_B$ from $z_S, \gamma_S$, $(z_{v,j})_{j=1,...,N_v}$ such that the normal velocity $\mathbf{u} \cdot \mathbf{n} = -\mathrm{Im}(\hat{\mathbf{u}}^{z_{S,e}}|_{z_{S,e}})$ satisfies the impermeability boundary condition on the bottom. This problem is then simpler than the previous one:

$$ \mathrm{pv} \int_{0}^{L_B} \gamma_B(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) \right] de' = \mathrm{RHS}_{VB}(e) \quad (3.1) $$

where

$$ \begin{aligned} \mathrm{RHS}_{VB}(e) &= - \int_{0}^{L_S} \gamma_S(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\ &\quad - \int_{\mathcal{D}_F} \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - y}{L/\pi} \right) \right] \omega_0 \, dy \\ &= - \int_{0}^{L_S} \gamma_S(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\ &\quad - \sum_{j=1}^{N_v} \gamma_{v,j} \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_{v,j}}{L/\pi} \right) \right] \\ &\quad - \frac{\omega_0}{4\pi} \int_{0}^{L_S} \ln \left( \cosh \mathrm{Im} \frac{z_B(e) - z_S(e')}{L/(2\pi)} - \cos \mathrm{Re} \frac{z_B(e) - z_S(e')}{L/(2\pi)} \right) \mathrm{Im} \left[ z_{B,e}(e) \overline{z_{S,e}(e')} \right] de' \\ &\quad + \frac{\omega_0}{4\pi} \int_{0}^{L_B} \ln \left( \cosh \mathrm{Im} \frac{z_B(e) - z_B(e')}{L/(2\pi)} - \cos \mathrm{Re} \frac{z_B(e) - z_B(e')}{L/(2\pi)} \right) \mathrm{Im} \left[ z_{B,e}(e) \overline{z_{B,e}(e')} \right] de', \end{aligned} $$

together with the circulation assumptions:

$$ \int_{0}^{L_B} \gamma_B(e') de' = -\gamma. $$

This problem is related to the vortex method in the case of an impermeable boundary. The invertibility of this problem was studied in details in [4], and the discretization will be described in Section 4.3.

This justifies Remark 1.2 in the case of the vortex formulation.
3.1.3. Initial $\mu_{S,0}$ for the dipole formulation. Regarding the dipole formulation, we will often have to construct $\mu_{S,0}$ knowing $g = \mathbf{u}_F \cdot \mathbf{n}|_{\Gamma_S}$. As usual, $z_B, z_S, \gamma, \omega_0$ and $(\gamma_{v,j}, z_{v,j})_{j=1,...,N_v}$ are initially given.

The first step is to construct $\mathbf{u}_{\omega,\gamma}$ if $(\omega, \gamma) \neq (0,0)$ (see Section 2.3). Given an expression of $\widehat{\mathbf{u}}_{\omega,\gamma}$, we are looking first for $\tilde{\gamma}_S, \tilde{\gamma}_B$ such that the associated $\mathbf{u}_R$ (2.18) solves the elliptic problem coming from (2.6), (2.7) and (2.10), i.e. we consider the unique solution of

$$ \begin{aligned} & \mathrm{pv} \int_0^{L_S} \tilde{\gamma}_S(e') \mathrm{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{d}e' \\ & + \int_0^{L_B} \tilde{\gamma}_B(e') \mathrm{Im} \left[ \frac{z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \right] \mathrm{d}e' = -g|z_S(e)| - \mathrm{Im} \left[ z_{S,e}(e) \widehat{\mathbf{u}}_{\omega,\gamma}(z_S(e)) \right] \end{aligned} $$

and

$$ \begin{aligned} & \int_0^{L_S} \tilde{\gamma}_S(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{d}e' \\ & + \mathrm{pv} \int_0^{L_B} \tilde{\gamma}_B(e') \mathrm{Im} \left[ \frac{z_{B,e}(e)}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) \right] \mathrm{d}e' = 0 \end{aligned} $$

together with the circulation assumptions

$$ \int_0^{L_S} \tilde{\gamma}_S(e') \mathrm{d}e' = \int_0^{L_B} \tilde{\gamma}_B(e') \mathrm{d}e' = 0. $$

Of course, if we have chosen the third construction of $\widehat{\mathbf{u}}_{\omega,\gamma}$, then we can replace $\mathrm{Im} [z_{S,e}(e)\widehat{\mathbf{u}}_{\omega,\gamma}(z_S(e))]$ by zero.

Finally, we set the initial value of $\mu_S$ as the anti-derivative of $\tilde{\gamma}_S$ with zero mean value

$$ \mu_{S,0}(e) = \int_0^e \tilde{\gamma}_S(e') \mathrm{d}e' - \frac{1}{L_S} \int_0^{L_S} \int_0^e \tilde{\gamma}_S(e') \mathrm{d}e' \mathrm{d}e. $$

3.1.4. Time dependent $\mu_B$ for the dipole formulation. When the time evolves, we need, for the dipole formulation, to construct $\mu_B$ from $z_S, \mu_S$, again from the initial data $z_B$. This problem is very simple, as we know that the potential obtained from $\mu_B$ and $\mu_S$ (see (2.20)) satisfies $\phi_B = 0$ in $D_B$. In particular, the limit of the potential (see Appendix A) by below $\Gamma_B$ vanishes, which reads

$$ \begin{aligned} & \frac{1}{2}\mu_B(e) + \int_0^{L_B} \mu_B(e') \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) z_{B,e}(e') \right] \mathrm{d}e' \\ & = - \int_0^{L_S} \mu_S(e') \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) z_{S,e}(e') \right] \mathrm{d}e' \quad (3.2) \end{aligned} $$
where the function in the left hand side integral is extended for $e = e'$ by

$$-\mu_B(e) \operatorname{Re} \left[ \frac{1}{4\pi i} \frac{z_{B,ee}(e)}{z_{B,e}(e)} \right].$$

By uniqueness of $\mu_B$ satisfying such an equation, we state that it is enough to solve it for $\mu_S$ given. This operator is different than the one in (3.1) and corresponds to the operator $A^*$ in [4]. We will highlight in Section 4.3 that it is possible to interpret this problem as a small perturbation of the identity, and its inverse will be obtained in the form of a Neumann series.

This justifies Remark 1.2 in the case of the dipole formulation.

*Remark 3.1.* This problem is easily stated as we are looking for $\mu_B$ such that $\phi_B = 0$ below the bottom, which is possible only if we have constructed $\mathbf{u}_{\omega,\gamma}$ tangent to the bottom (2.10).

## 3.2. Displacement of the free surface.
For the vortex formulation, we already have $\gamma_S$ and $\gamma_B$, whereas for the dipole formulation we simply construct them from $\mu_S$ and $\mu_B$: $\tilde{\gamma}_B = \partial_e \mu_B(e)$ and $\tilde{\gamma}_S = \partial_e \mu_S(e)$.

With $\gamma_S$ and $\gamma_B$, it is now possible from (2.14) to compute the velocity anywhere. We recall that the tangential velocity is discontinuous at the water-air interface. We can thus assume that the free surface is moving as

$$\partial_t \mathbf{z}_S(t, e) = \alpha \mathbf{u}_F(t, z_F(e)) + (1 - \alpha) \mathbf{u}_A(t, z_F(e))$$

with a parameter $\alpha \in [0, 1]$. By continuity of the normal component, the evolution of the free surface does not depend on the choice of $\alpha$. It is however an interesting parameter from a numerical point of view. It controls the distribution of points on the interface. A choice of $\alpha$ can, for example, allow us to vary the resolution in space. In practice, we found that $\alpha = 1$ efficiently concentrates computational points near the tip of a breaking wave.

The limit formulae of Appendix A give

$$\begin{aligned}
\overline{\partial_t z_S(t, e)} &= \int_0^{L_S} \frac{\gamma_S(e') z_{S,e}(e) - \gamma_S(e) z_{S,e}(e')}{z_{S,e}(e)} \frac{1}{2Li} \cot\left(\frac{z_S(e) - z_S(e')}{L/\pi}\right) de' \\
&\quad + \int_0^{L_B} \gamma_B(e') \frac{1}{2Li} \cot\left(\frac{z_S(e) - z_B(e')}{L/\pi}\right) de' \\
&\quad + \frac{2\alpha - 1}{2} \frac{\gamma_S(e)}{z_{S,e}(e)} + \sum_{j=1}^{N_v} \frac{\gamma_{v,j}}{2Li} \cot\left(\frac{z_S(e) - z_{v,j}}{L/\pi}\right) \\
&\quad + \frac{\omega_0}{4\pi} \int_0^{L_S} \ln\left(\cosh \operatorname{Im} \frac{z_S(e) - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_S(e')}{L/(2\pi)}\right) \overline{z_{S,e}(e')} de' \\
&\quad - \frac{\omega_0}{4\pi} \int_0^{L_B} \ln\left(\cosh \operatorname{Im} \frac{z_S(e) - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_B(e')}{L/(2\pi)}\right) \overline{z_{B,e}(e')} de'.
\end{aligned}
\quad (3.3)$$
This formula highlights that the dependence on $\alpha$ appears only through a tangent vector field $\frac{\gamma_S(e)}{|z_{S,e}(e)|}z_{S,e}(e)$.

The displacement of the point vortices is an obvious application of this formula. For all $i = 1, \dots, N_v$,

$$ \partial_t \overline{z_{v,i}(t)} = \int_0^{L_S} \gamma_S(e') \frac{1}{2Li} \cot \left( \frac{z_{v,i} - z_S(e')}{L/\pi} \right) de' + \int_0^{L_B} \gamma_B(e') \frac{1}{2Li} \cot \left( \frac{z_{v,i} - z_B(e')}{L/\pi} \right) de' \\ + \sum_{j=1, j \neq i}^{N_v} \frac{\gamma_{v,j}}{2Li} \cot \left( \frac{z_{v,i} - z_{v,j}}{L/\pi} \right) \\ + \frac{\omega_0}{4\pi} \int_0^{L_S} \ln \left( \cosh \operatorname{Im} \frac{z_{v,i} - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_{v,i} - z_S(e')}{L/(2\pi)} \right) \overline{z_{S,e}(e')} de' \\ - \frac{\omega_0}{4\pi} \int_0^{L_B} \ln \left( \cosh \operatorname{Im} \frac{z_{v,i} - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_{v,i} - z_B(e')}{L/(2\pi)} \right) \overline{z_{B,e}(e')} de'. $$

In the case of the dipole formulation, we simply use (2.18) to get $\mathbf{u}_R$ everywhere. We thus need the expression of $\mathbf{u}_{\omega,\gamma}$, see Section 3.1.3 for these formulae depending on the physical setting. In particular, we notice that often, we do not have the velocity in the air, hence we should only consider the case

$$ \partial_t \mathbf{z}_S(t, e) = \mathbf{u}_F(t, z_F(e)) $$

which makes sense for the single-fluid water-waves equations. We thus write

$$ \partial_t \overline{z_S(t, e)} = \int_0^{L_S} \frac{\tilde{\gamma}_S(e') z_{S,e}(e) - \tilde{\gamma}_S(e) z_{S,e}(e')}{z_{S,e}(e)} \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) de' \\ + \int_0^{L_B} \tilde{\gamma}_B(e') \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) de' + \frac{1}{2} \frac{\tilde{\gamma}_S(e)}{z_{S,e}(e)} + \widehat{\mathbf{u}}_{\omega,\gamma}(z_S(e)), $$

where the last formula of $\widehat{\mathbf{u}}_{\omega,\gamma}$ in Section 3.1.3 has to be desingularized when $x = z_S(e)$ as we did in the previous formula (3.3).

It is worth stressing that some freedom is left on the choice of $\alpha$, which will not affect the shape of the solution, but only the tangential distribution of points on the interface.

In the case of the bi-fluid problem, we have constructed $\mathbf{u}_{\omega,\gamma}$ tangent to free surface, it is thus more natural to include the parameter $\alpha$

$$ \partial_t \overline{z_S(t, e)} = \int_0^{L_S} \frac{\tilde{\gamma}_S(e') z_{S,e}(e) - \tilde{\gamma}_S(e) z_{S,e}(e')}{z_{S,e}(e)} \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) de' + \int_0^{L_B} \tilde{\gamma}_B(e') \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) de' + \frac{2\alpha - 1}{2} \frac{\tilde{\gamma}_S(e)}{z_{S,e}(e)} + \alpha \widehat{\mathbf{u}}_{\omega,\gamma}(z_S(e)). \quad (3.4) $$
3.3. **Bernoulli equation and dipole formulation.** Writing $\mathbf{u} = \mathbf{u}_{\omega,\gamma} + \nabla\phi_F$, the Euler equations take the form

$$ \nabla\left[\partial_t\phi_F + \frac{1}{2}|\mathbf{u}|^2\right] + \partial_t\mathbf{u}_{\omega,\gamma} + (\mathrm{curl}\,\mathbf{u})\mathbf{u}^{\perp} = -\nabla\left[\frac{p_F}{\rho_F} + gx_2\right] $$

where $\rho_F$ is the density of the fluid and $g$ is the gravity acceleration. Hence in the neighborhood of the free surface where the vorticity is constant, the Euler equations can be reduced to the modified Bernoulli equation

$$ \partial_t\phi_F + \frac{1}{2}|\nabla\phi_F + \mathbf{u}_{\omega,\gamma}|^2 + \phi_{t,\omega,\gamma} - \omega_0(\psi_{\omega,\gamma} + \tilde{\psi}_F) = -\frac{p_F}{\rho_F} - gx_2, $$

with $\psi_{\omega,\gamma}$ the stream function of $\mathbf{u}_{\omega,\gamma}$ and $\phi_{t,\omega,\gamma}$ the potential of $\partial_t\mathbf{u}_{\omega,\gamma}^4$.

*Remark 3.2.* For all the examples given above, the difficult component to express is the stream function associated to the constant part. This is because it is not obvious how to express $\int_{\mathcal{D}_F} G(\mathbf{x}-\mathbf{y}) \, dy$ as an integral over the boundaries. Hence, it is easier to assume $\omega_0 = 0$ which removes the presence of $\psi_{\omega,\gamma}$.

We should also observe that the potential $\phi_{t,\omega,\gamma}$ is even more complicated to obtain. Therefore, in the sequel, we only consider stationary $\mathbf{u}_{\omega,\gamma}$ where we can forget $\partial_t\mathbf{u}_{\omega,\gamma}$ and $\phi_{t,\omega,\gamma}$. Of course, this implies that we also assume $\gamma_{v,j} = 0$, hence $\omega = 0$ and we can replace $\mathbf{u}_{\omega,\gamma}$ with $\mathbf{u}_{\gamma}$.

From now on, we will restrict our attention to the following two situations for the dipole formulation:

* the bi-fluid and the single fluid water-waves equations without circulation and vorticity, (i.e. $\mathbf{u}_{\gamma} = 0$);

* the single fluid water-waves equations with circulation and without vorticity, (i.e. $\mathbf{u}_{\gamma} = \gamma\mathbf{e}_1/L$ for the flat bottom and $\mathbf{u}_{\gamma} = \mathbf{H}$, initially constructed for other bottom).

For the bi-fluid water-waves equations, we also consider the Bernoulli equation in the air

$$ \partial_t\phi_A + \frac{1}{2}|\nabla\phi_A|^2 = -\frac{p_A}{\rho_A} - gx_2. $$

For the single-fluid water-waves equations, the density of the air is neglected and the pressure in the air is constant.

It is useful to count here boundary conditions for the fluid domain $\mathcal{D}_F$. At the bottom boundary $\Gamma_B$ the normal component of the flow vanishes, which provides the needed boundary condition. At the top boundary $\Gamma_S$, in the bi-fluid problem both the normal component of velocity and the pressure are continuous with that is the air domain $\mathcal{D}_A$ above, across the boundary $\Gamma_S$. These two continuity relations are thus enough to close the fluid system in $\mathcal{D}_F$ (morally the number of jump conditions needed at the boundary between two domains is $n_1 + n_2$, where $n_1$ and $n_2$ are the numbers of outer conditions

<footnote n="4">⁴Such a potential exists in the neighborhood of the free surface, because $\mathrm{curl}\,\partial_t\mathbf{u}_{\omega,\gamma} = \partial_t\omega_0 = 0$ and by the conservation of circulation.</footnote>
required for the PDE solution in domains 1 and 2 respectively). The single-fluid water-waves problem, corresponds to the limit of a vanishing density for the air. The two jump conditions on $\Gamma_S$ are then replaced by a single boundary condition on pressure $p = 0$, which again is enough to close the Euler system in $\mathcal{D}_F$.

In any case, we need to relate $p_A$ and $p_F$. In the presence of surface tension, the pressure is not continuous and a pressure jump is achieved, which is directly related to the surface curvature $\kappa$

$$[p\mathbf{n}] = (p_A - p_F)\mathbf{n} = -\sigma\kappa\mathbf{n}. \qquad (3.5)$$

The surface tension coefficient $\sigma$ is here related to capillary effects (e.g., [27, §9.1.2]).

Setting the atmospheric pressure to 0 for the single-fluid equations, we consider the limit of the Bernoulli equation at the interface

$$\begin{aligned} \partial_t (\phi_F(z_S(e))) - \partial_t \mathbf{z}_S(e) \cdot (\nabla \phi_F)(z_S(e)) + \frac{1}{2} |\nabla \phi_F + \mathbf{u}_\gamma|^2(z_S(e)) \\ = -\frac{\sigma}{\rho_F} \kappa(z_S(e)) - g z_{S,2}(e), \end{aligned}$$

hence we get

$$\begin{aligned} \frac{1}{2} \partial_t \mu_S(e) &= \frac{1}{2} \partial_t (\phi_F(z_S(e))) - \frac{1}{2} \partial_t (\phi_A(z_S(e))) = \partial_t (\phi_F(z_S(e))) - \frac{1}{2} \partial_t \Phi_S(e) \\ &= -\frac{1}{2} \partial_t \Phi_S(e) + \partial_t \mathbf{z}_S(e) \cdot (\nabla \phi_F)(z_S(e)) - \frac{1}{2} |\nabla \phi_F + \mathbf{u}_\gamma|^2(z_S(e)) \quad (3.6) \\ &\quad - \frac{\sigma}{\rho_F} \kappa(z_S(e)) - g z_{S,2}(e). \end{aligned}$$

For the bi-fluid water-waves equation without circulation, we need to consider the Bernoulli equation in the air, hence performing the difference we get

$$\begin{aligned} \partial_t \mu_S(e) + \partial_t \mathbf{z}_S(e) \cdot (\nabla \phi_A - \nabla \phi_F)(z_S(e)) + \frac{1}{2} (|\nabla \phi_F|^2 - |\nabla \phi_A|^2)(z_S(e)) \\ = \frac{\rho_F - \rho_A}{\rho_F \rho_A} p_F(z_S(e)) - \frac{\sigma}{\rho_A} \kappa(z_S(e)), \end{aligned}$$

whereas the sum gives

$$\begin{aligned} \partial_t \Phi_S(e) - \partial_t \mathbf{z}_S(e) \cdot (\nabla \phi_A + \nabla \phi_F)(z_S(e)) + \frac{1}{2} (|\nabla \phi_F|^2 + |\nabla \phi_A|^2)(z_S(e)) \\ = -\frac{\rho_F + \rho_A}{\rho_F \rho_A} p_F(z_S(e)) + \frac{\sigma}{\rho_A} \kappa(z_S(e)) - 2g z_{S,2}(e). \end{aligned}$$

We remove the pressure by multiplying the second equation by the Atwood number

$$A_{tw} = \frac{\rho_F - \rho_A}{\rho_F + \rho_A} \qquad (3.7)$$

and finally obtain

$$\begin{aligned} \frac{1}{2} \partial_t \mu_S(e) = & -\frac{A_{tw}}{2} \partial_t \Phi_S(e) + \frac{1}{2} \partial_t \mathbf{z}_S(e) \cdot ((A_{tw} + 1)\nabla \phi_F + (A_{tw} - 1)\nabla \phi_A)(z_S(e)) \\ & -\frac{1}{4} \left( (A_{tw} + 1)|\nabla \phi_F|^2 + (A_{tw} - 1)|\nabla \phi_A|^2 \right) (z_S(e)) \end{aligned}$$
$$ -\frac{(A_{tw} + 1)\sigma}{2\rho_F} \kappa(z_S(e)) - gA_{tw}z_{S,2}(e) \qquad (3.8) $$

because $\frac{A_{tw}-1}{\rho_A} = \frac{-2}{\rho_F+\rho_A} = -\frac{A_{tw}+1}{\rho_F}$. Even if the derivation differs, it is worth stressing that the single-fluid water-waves equations without circulation are recovered when setting $A_{tw} = 1$ in the above equation.

We have already derived the expression for $\partial_t z_S(e)$ and in the similar way for $\overline{\nabla\phi_F}(z_S(e))$ and $\overline{\nabla\phi_A}(z_S(e))$, our next step is to compute $\partial_t \Phi_S(e)$. From (2.23)

$$ \frac{\partial_t \Phi_S(e)}{2} = \int_0^{L_S} (\partial_t \mu_S(e') - \partial_t \mu_S(e)) \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) z_{S,e}(e') \right] \mathrm{d}e' \\ + \int_0^{L_B} \partial_t \mu_B(e') \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) z_{B,e}(e') \right] \mathrm{d}e' \\ + \int_0^{L_S} (\mu_S(e) - \mu_S(e')) \mathrm{Re} \left[ \frac{\pi}{2L^2i} \sin^{-2} \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) (\partial_t z_S(e) - \partial_t z_S(e')) z_{S,e}(e') \right] \mathrm{d}e' \\ + \int_0^{L_S} (\mu_S(e') - \mu_S(e)) \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \partial_t z_{S,e}(e') \right] \mathrm{d}e' \\ - \int_0^{L_B} \mu_B(e') \mathrm{Re} \left[ \frac{\pi}{2L^2i} \sin^{-2} \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \partial_t z_S(e) z_{B,e}(e') \right] \mathrm{d}e' \qquad (3.9) $$

where the third right hand side integral concerns a continuous function whose value for $e' = e$ is

$$ \mu_{S,e}(e) \mathrm{Re} \left[ \frac{1}{2\pi i} \frac{\partial_t z_{S,e}(e)}{z_{S,e}(e)} \right] $$

whereas in the fourth integral the extension is

$$ -\mu_{S,e}(e) \mathrm{Re} \left[ \frac{1}{2\pi i} \frac{\partial_t z_{S,e}(e)}{z_{S,e}(e)} \right] $$

which is exactly opposite to the first one, and then can be omitted.

Finally, substituting (3.9) into (3.6) or (3.8), yields an equation of the form

$$ A_S^*[\partial_t \mu_S](e) + C_D[\partial_t \mu_B](e) = G_{D,1}(e), $$

where the operator $A_S$ is the same kind of operator as in (3.2) (see Section 4.3 for an explanation on how to discretize such an operator, the discrete expressions of $A_S$, $C_D$ and $G_{D,1}$ are given in Appendix B).
As the above equation involves $\partial_t \mu_B$, we derive another equation differentiating (3.2) with respect to time:

$$ 
\begin{aligned}
& \int_0^{L_S} \partial_t \mu_S(e') \operatorname{Re} \left[ \frac{z_{S,e}(e')}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\
& + \frac{1}{2} \partial_t \mu_B(e) + \int_0^{L_B} \partial_t \mu_B(e') \operatorname{Re} \left[ \frac{z_{B,e}(e')}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) \right] de' \\
& = - \int_0^{L_S} \mu_S(e') \operatorname{Re} \left[ \frac{\partial_t z_{S,e}(e')}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\
& - \int_0^{L_S} \mu_S(e') \operatorname{Re} \left[ \frac{\pi z_{S,e}(e') \partial_t z_S(e')}{2L^2i} \sin^{-2} \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de'
\end{aligned}
$$

which gives an equation of the form

$$ D_D[\partial_t \mu_S](e) + A_B^*[\partial_t \mu_B](e) = G_{D,2}(e) $$

where $A_B^*$ is precisely the same operator as on the left hand side of (3.2), which was already inverted. From this equation, we have

$$ \partial_t \mu_B = A_B^{*-1} [G_{D,2} - D_D[\partial_t \mu_S]] $$

which means that $\partial_t \mu_S$ will be obtained by solving

$$ (A_S^* - C_D A_B^{*-1} D_D)[\partial_t \mu_S] = G_{D,1} - C_D A_B^{*-1} G_{D,2}. \quad (3.10) $$

We will discuss in Section 4.3 that $A_S^* - C_D A_B^{*-1} D_D$ can be seen as a perturbation of a simple matrix, though not the identity. The discrete version is also given in Appendix B.

This equation, together with (3.4) corresponds to the vortex formulation and the first part of Theorem 1.3.

### 3.4. Euler equation and vortex formulation.
In [6], the authors differentiate the equation for $\partial_t \mu$ (3.8) to get the equation for $\partial_t \gamma$. This is difficult to justify because the kernels in the integrals are singular. Let us also stress that such a derivation yields a vortex formation which should only be used without circulation, i.e. with zero mean value for $\gamma_{S,0}$.

Alternatively, we write the Euler equations in $\mathcal{D}_F$

$$ \partial_t \mathbf{u}_F + (\mathbf{u}_F \cdot \nabla) \mathbf{u}_F = \frac{-\nabla p_F}{\rho_F} - g\mathbf{e}_2. $$

For the single-fluid water-waves equation, we have $\nabla p_A = 0$, whereas for the bi-fluid water-waves equations we also write the Euler equations in $\mathcal{D}_A$

$$ \partial_t \mathbf{u}_A + (\mathbf{u}_A \cdot \nabla) \mathbf{u}_A = \frac{-\nabla p_A}{\rho_A} - g\mathbf{e}_2. $$
As for the dipole formulation, we need to relate the pressures on both sides of the interface using the continuity of the normal component of the stress tensor at the interface (3.5). This implies by differentiating with respect to e

$$ \mathbf{z}_{S,e}(e) \cdot (\nabla p_A - \nabla p_F)(z_S(e)) = -\sigma \frac{d}{de} \left( \kappa(z_S(e)) \right). $$

So we need to consider the limit of the tangential part of the Euler equations at the interface.

For the single-fluid water-waves equations, one can simply replace $\mathbf{z}_{S,e}(e) \cdot \nabla p_F(z_S(e))$ by $\sigma \frac{d}{de} \left( \kappa(z_S(e)) \right)$ and obtain

$$ \partial_t (\mathbf{u}_F(z_S(e)) \cdot \mathbf{z}_{S,e}(e)) + \left[ \left( (\mathbf{u}_F(z_S(e)) - \partial_t \mathbf{z}_S(e)) \cdot \nabla \right) \mathbf{u}_F(z_S(e)) \right] \cdot \mathbf{z}_{S,e}(e) \\ - \mathbf{u}_F(z_S(e)) \cdot \partial_t \mathbf{z}_{S,e}(e) = - \frac{\sigma}{\rho_F} \frac{d}{de} \left( \kappa(z_S(e)) \right) - g z_{S,e,2}. $$

In order to introduce $\partial_t \gamma_S$, we use

$$ \begin{aligned} \Psi_S(e) &:= (\mathbf{u}_F + \mathbf{u}_A)(z_S(e)) \cdot \mathbf{z}_{S,e}(e) \\ &= \text{pv} \int_0^{L_S} \gamma_S(e') \text{Re} \left[ \frac{z_{S,e}(e)}{Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{d}e' \\ &\quad + \int_0^{L_B} \gamma_B(e') \text{Re} \left[ \frac{z_{S,e}(e)}{Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \right] \mathrm{d}e' \\ &\quad + \sum_{j=1}^{N_v} \gamma_{v,j} \text{Re} \left[ \frac{z_{S,e}(e)}{Li} \cot \left( \frac{z_S(e) - z_{v,j}}{L/\pi} \right) \right] \\ &\quad + \frac{\omega_0}{2\pi} \int_0^{L_S} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_S(e')}{L/(2\pi)} \right) \operatorname{Re} \left[ z_{S,e}(e) \overline{z_{S,e}(e')} \right] \mathrm{d}e' \\ &\quad - \frac{\omega_0}{2\pi} \int_0^{L_B} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_B(e')}{L/(2\pi)} \right) \operatorname{Re} \left[ z_{S,e}(e) \overline{z_{B,e}(e')} \right] \mathrm{d}e' \end{aligned} $$

so by the definition of $\gamma_S$ (2.12), we write

$$ \begin{aligned} \frac{1}{2} \partial_t \gamma_S(e) &= \partial_t (\mathbf{u}_F(z_S(e)) \cdot \mathbf{z}_{S,e}(e)) - \frac{1}{2} \partial_t \Psi_S(e) \\ &= -\frac{1}{2} \partial_t \Psi_S(e) - \left[ \left( (\mathbf{u}_F(z_S(e)) - \partial_t \mathbf{z}_S(e)) \cdot \nabla \right) \mathbf{u}_F(z_S(e)) \right] \cdot \mathbf{z}_{S,e}(e) \quad (3.11) \\ &\quad + \mathbf{u}_F(z_S(e)) \cdot \partial_t \mathbf{z}_{S,e}(e) - \frac{\sigma}{\rho_F} \frac{d}{de} \left( \kappa(z_S(e)) \right) - g z_{S,e,2}. \end{aligned} $$

For the bi-fluid formulation, we proceed as for the dipole formulation, i.e.

• we compute from both Euler equations

$$ \partial_t \gamma_S(e) = \partial_t (\mathbf{u}_F(z_S(e)) \cdot \mathbf{z}_{S,e}(e)) - \partial_t (\mathbf{u}_A(z_S(e)) \cdot \mathbf{z}_{S,e}(e)); $$
* we express in the same way

$$ \partial_t \Psi_S(e) = \partial_t (\mathbf{u}_F(z_S(e)) \cdot \mathbf{z}_{S,e}(e)) + \partial_t (\mathbf{u}_A(z_S(e)) \cdot \mathbf{z}_{S,e}(e)); $$

* we replace $\mathbf{z}_{S,e}(e) \cdot \nabla p_A(z_S(e))$ by $\mathbf{z}_{S,e}(e) \cdot \nabla p_F(z_S(e)) - \sigma \frac{d}{de}(\kappa(z_S(e)))$;

* we remove the pressure term by multiplying the second equation by the Atwood number (3.7), and by adding the two equations, and we use $\frac{A_{tw}-1}{\rho_A} = \frac{-2}{\rho_F+\rho_A} = -\frac{A_{tw}+1}{\rho_F}$.

In the end, we get a slightly modified equation for $\partial_t \gamma_S(e)$:

$$ \begin{aligned} \frac{1}{2} \partial_t \gamma_S(e) = & -\frac{A_{tw}}{2} \partial_t \Psi_S(e) - \frac{1+A_{tw}}{2} \left[ \left( (\mathbf{u}_F(z_S(e)) - \partial_t \mathbf{z}_S(e)) \cdot \nabla \right) \mathbf{u}_F(z_S(e)) \right] \cdot \mathbf{z}_{S,e}(e) \\ & + \frac{1-A_{tw}}{2} \left[ \left( (\mathbf{u}_A(z_S(e)) - \partial_t \mathbf{z}_S(e)) \cdot \nabla \right) \mathbf{u}_A(z_S(e)) \right] \cdot \mathbf{z}_{S,e}(e) - g A_{tw} z_{S,e,2} \\ & + \left( \frac{1+A_{tw}}{2} \mathbf{u}_F - \frac{1-A_{tw}}{2} \mathbf{u}_A \right) (z_S(e)) \cdot \partial_t \mathbf{z}_{S,e}(e) - \frac{(1+A_{tw})\sigma}{2\rho_F} \frac{d}{de}(\kappa(z_S(e))) \end{aligned} \quad (3.12) $$

which coincides with the first one when we set $A_{tw} = 1$.

We may simplify the second and third right hand side terms by

$$ \begin{aligned} \left( (\mathbf{u}_F(z_S) - \partial_t \mathbf{z}_S) \cdot \nabla \right) \mathbf{u}_F(z_S) &= (1 - \alpha) \frac{\gamma_S}{|z_{S,e}|^2} \left( \mathbf{z}_{S,e} \cdot \nabla \right) \mathbf{u}_F(z_S) \\ &= (1 - \alpha) \frac{\gamma_S}{|z_{S,e}|^2} \partial_e \left( \mathbf{u}_F(z_S) \right), \\ \left( (\mathbf{u}_A(z_S) - \partial_t \mathbf{z}_S) \cdot \nabla \right) \mathbf{u}_A(z_S) &= -\alpha \frac{\gamma_S}{|z_{S,e}|^2} \left( \mathbf{z}_{S,e} \cdot \nabla \right) \mathbf{u}_A(z_S) = -\alpha \frac{\gamma_S}{|z_{S,e}|^2} \partial_e \left( \mathbf{u}_A(z_S) \right). \end{aligned} $$

We then substitute the computation of $\partial_t \Psi_S$ (see Appendix C) into (3.11) or (3.12) to get an equation of the form

$$ A_S[\partial_t \gamma_S](e) + C_V[\partial_t \gamma_B](e) = G_{V,1}(e) $$

where $A_S$ is the adjoint operator of $A_S^*$ in the dipole formulation, see Section 4.3 for properties of these operators.

Finally, we differentiate with respect to time (3.1) to get an equation of the form

$$ D_V[\partial_t \gamma_S](e) + B_B[\partial_t \gamma_B](e) = G_{V,2}(e) $$

where the conservation of circulation $\int \partial_t \gamma_B = 0$ is included in the last line of the discretization. This allows to obtain $\partial_t \gamma_S$ by solving

$$ (A_S - C_V B_B^{-1} D_V)[\partial_t \gamma_S] = G_{V,1} - C_V B_B^{-1} G_{V,2}. \quad (3.13) $$

Again, the discrete forms are given in Appendix C.

This equation, together with (3.3) corresponds to the vortex formulation and the last part of Theorem 1.3.
Remark 3.3. In [4], we have developed a method referred to as the “fluid charge method”. In this method, after retrieving $\mathbf{u}_{\omega,\gamma}$ we have written $\mathbf{u}_R = \nabla\phi_F$ and solved the Laplace problem with homogeneous Neumann boundary condition on $\Gamma_B$. This method relies on an extension of $\phi$ in $\mathcal{D}_B$ by continuity. We can establish that

$$ \widehat{\nabla}\phi(x) = \int_{0}^{L_S} \sigma_S(e) \frac{1}{2L} \cot\left(\frac{x - z_S(e)}{L/\pi}\right) \mathrm{d}e + \int_{0}^{L_B} \sigma_B(e) \frac{1}{2L} \cot\left(\frac{x - z_B(e)}{L/\pi}\right) \mathrm{d}e, $$

where

$$ \sigma_S(e) := |z_{S,e}(e)| \left[ \lim_{z \in \mathcal{D}_A \to z_S(e)} \partial_n \phi - \lim_{z \in \mathcal{D}_F \to z_S(e)} \partial_n \phi \right] $$

$$ \sigma_B(e) := |z_{B,e}(e)| \left[ \lim_{z \in \mathcal{D}_F \to z_B(e)} \partial_n \phi - \lim_{z \in \mathcal{D}_B \to z_B(e)} \partial_n \phi \right] = -|z_{B,e}(e)| \left[ \lim_{z \in \mathcal{D}_B \to z_B(e)} \partial_n \phi \right]. $$

Indeed, in this formulation, the tangential part is continuous whereas the normal part has a jump. Therefore, we can adapt Section 3.1 to find $\sigma_B$ such that $\mathbf{u}_R$ satisfies the impermeability condition on $\Gamma_B$. Note that this problem is related to the inversion of the operator $A$ in [4]. Next, we can use the previous formula to get the displacement of the free surface in terms of $\sigma_S$.

Unfortunately, in the case of water-waves, we do not have an equation for $\partial_t \sigma_S$. In Sections 3.3 and 3.4, we have used the continuity of the normal component of the stress tensor (3.5). For the dipole formulation, this equation on $p_F - p_A$ provides the connexion between the two Bernoulli equations and allows to get an equation on $\partial_t \mu_S = \frac{1}{2}\partial_t(\phi_F(z_S(e))) - \frac{1}{2}\partial_t(\phi_A(z_S(e)))$. For the vortex formulation, a differentiation along the free surface of this relation on $p_F - p_A$ yields an equation on $\mathbf{z}_{S,e}(e) \cdot (\nabla p_A - \nabla p_F)(z_S(e))$. This establishes a connexion between the tangential part of the two Euler equations. Here, it is important that $\gamma_S$ corresponds to the jump of the tangential velocities. In the fluid charge method, $\sigma_S$ corresponds to the jump of the normal velocities. To get an equation for $\partial_t \sigma_S$, we would thus need a relation for $\mathbf{z}_{S,e}(e)^\perp \cdot (\nabla p_A - \nabla p_F)(z_S(e))$, i.e. a sort of continuity of the normal derivative of the normal component of the stress tensor, which has no physical meaning. Note that, defining the velocity above the free surface such that the normal component has a jump implies that it does not corresponds to the air velocity.

Let us note that Baker in [6, Equation (4.4)] considered either the vortex or dipole formulation for the free surface combined with the fluid charge method at the bottom. Because of the extension of $\psi$ above the fluid and $\phi$ below the fluid domain, Proposition 2.1 cannot be applied to such a formulation. For this reason, we have chosen in the previous sections to consider the same formulation for the free surface and for the bottom.

### 3.5. The deep-water case.
It is easy to derive a deep-water formulation removing contributions from the bottom. Following line by line the previous sections without the presence of $\mathcal{D}_B$ and $\Gamma_B$, i.e. assuming that the fluid domain is infinite in the vertical direction, we can get the following model

* the dipole formulation for the bi-fluid water-waves equations without vorticity and circulation, where the velocity is given by $\gamma_S = \partial_t \mu_S$. The equation verified by $\partial_t \mu_S$ stays exactly (3.8), but dropping all terms involving $\mu_B$ in the expression
(3.9) for $\partial_t\Phi_S$. This yields

$$A_S^*[\partial_t\mu_S](e) = G_{D,1}(e),$$

where $G_{D,1}$ is giving in (B.1) where we remove $\mathbf{u}_\gamma$ and $\mu_B$;

* the dipole formulation for the single-fluid water-waves equations with circulation and without vorticity, where $\mathbf{u}_\gamma = \frac{\gamma}{L}\mathbf{e}_1$ and $\tilde{\gamma}_S = \partial_e\mu_S$. The equation verified by $\partial_t\mu_S$ stays exactly (3.6), but again dropping all terms involving $\mu_B$ in the expression (3.9) for $\partial_t\Phi_S$. This yields

$$A_S^*[\partial_t\mu_S](e) = G_{D,1}(e),$$

where $G_{D,1}$ is giving in (B.1) where we remove $\mu_B$ and replace $A_{tw}$ by 1 and $\mathbf{u}_\gamma$ by $\frac{\gamma}{L}\mathbf{e}_1$;

* the vortex formulation with circulation and where the vorticity is composed of point vortices (no constant part). The velocity is given by $\gamma_S$ and the equation verified by $\partial_t\gamma_S$ stays exactly (3.12), but dropping all terms involving $\gamma_B$ in the expression for $\partial_t\Psi_S$, see Appendix C. This yields

$$A_S[\partial_t\gamma_S](e) = G_{V,1}(e).$$

Let us note that the equations obtained in this case are very close to [6, Equations (2.14)-(2.17)], but where we have used the desingularization (2.23), which allows us to justify the derivatives and to handle only classical integrals. In this earlier article, principal value integrals with $x/\sin^2 x$ singularities may be a cause of numerical instabilities.

## 4. NUMERICAL DISCRETIZATION

### 4.1. Time integration.
Whereas most earlier numerical studies on water waves breaking used high-order Runge-Kutta integrators [3, 7, 33], we preferred to restrict our study to a second order in time, but symplectic integrator for harmonic oscillators. We use the so-called Verlet integrator, which amounts to using a staggered grid in time, and preserves the Hamiltonian structure of harmonic oscillators.

The governing equations take the form

$$\partial_t X = G(Z, X), \quad \partial_t Z = F(Z, X),$$

where $X \equiv \gamma$ in the vortex formulation and $X \equiv \mu$ in the dipole formulation. $F$ and $G$ denote here non-linear differential operators. These are discretized in the form

$$X^{n+1/2} - X^{n-1/2} = \Delta t G(Z^n, X^n), \quad Z^{n+1} - Z^n = \Delta t F(Z^{n+1/2}, X^{n+1/2}).$$

The right-hand-side of the first equation involves $X^n$ which is not known and that of the second equation similarly involves the unknown $Z^{n+1/2}$. These are respectively constructed as $2X^n \simeq X^{n+1/2} + X^{n-1/2}$ and $2Z^{n+1/2} \simeq Z^{n+1} + Z^n$ and calculated using a fixed point relaxation.

We observed numerically that this symplectic integrator offers better stability properties than standard Runge-Kutta integrator and yields remarkable conservation properties on test cases (see Section 5).
4.2. **Shifted grids in space.** Besides the use of a staggered mesh in time, a staggered mesh in space can also be used. This approach was for example used and fully justified (via a mathematical demonstration) in [4] to enforce impermeability boundary conditions.

Shifted grids are also used here for the free surface. Our aim is to avoid regularization techniques and yet desingularize the integrals involved in the computation. We reformulated all singular integrals as regular ones using relation (2.16) (see also Appendix A). The resulting integral is now non-singular, but can only be defined at the former singularity as a continuous prolongation. This extension necessarily involves higher order derivatives, which can induce some numerical errors when the curvature becomes large (i.e. in a situation relevant to wave breaking). Evaluating the integrals on a shifted dual grid resolves this problem as the function is at worst evaluated half a grid point away from the former singularity. The integral eventually needs to be interpolated on the original grid for time stepping. This introduces some numerical smoothing, which is however entirely controlled by the grid size and thus vanishes in the limit of a large number of points. Finally, all derivatives in space, with respect to $e$, in the discrete expressions are evaluated thanks to second order finite difference formula.

4.3. **Discretization and inversion of singular operator by Neumann series.** One of the main numerical difficulties lies in the resolution of linear systems with matrices related to singular kernel operators. If it is well known that the continuous operators are indeed invertible [16], the relevant discretization, the invertibility and the convergence of the discrete operators have recently been studied in [4]. The notation $A, B$ and their adjoints $A^*, B^*$ come from Equation (3.1) in this article, where the relation and the inversion is based on the Poincaré-Bertrand formula concerning the inversion of Cauchy integrals, see [4, Section 3] for full details.

Given an arc-length parametrization : $[0, L_B] \to \Gamma_B$, if

$$0 = e_{B,1} < e_{B,2} < \cdots < e_{B,N} < L_B$$

are close to the uniform distribution $\theta_i = (i - 1)L_B/N_B$, then the matrix $A_{B,N}^*$ appearing to compute $\mu_B$, see Section 3.1.4, defined as

$$A_{B,N}^*(i,j) = \frac{L_B}{N_B} \operatorname{Re} \left[ \frac{1}{2L_B i} \cot \left( \frac{z_B(e_{B,i}) - z_B(e_{B,j})}{L/\pi} \right) z_{B,e}(e_{B,j}) \right] \quad \forall i \neq j \in [1, N_B] \times [1, N_B],$$

$$A_{B,N}^*(i,i) = \frac{1}{2} - \frac{L_B}{N_B} \operatorname{Re} \left[ \frac{1}{2\pi i} \frac{z_{B,ee}(e_{B,i})}{z_{B,e}(e_{B,i})} \right] \quad \forall i \in [1, N_B],$$

is invertible and can be seen as a perturbation of $\frac{1}{2}I_2$. It is thus a well conditioned matrix, see for instance [4, Section 8.3]. It can be inverted very efficiently by a Neumann series. We refer to [4, Section 8], in particular Theorem 8.1 therein where the convergence rate is given.

Namely, we write

$$A_{B,N}^* = \frac{1}{2}(I_N - R_{B,N}), \quad \|R_{B,N}\| < 1,$$
which implies that

$$A_{B,N}^{*-1} = 2(I_N - R_{B,N})^{-1} = 2 \sum_{k=0}^{+\infty} R_{B,N}^k.$$

In view of a fixed point procedure, we denote $R_{B,N} := -2A_{B,N}^* + I_N$, and

$$U_{n+1} = 2 \sum_{k=0}^{n+1} R_{B,N}^k = R_{B,N} \left( 2 \sum_{k=0}^{n} R_{B,N}^k \right) + 2I_N = R_{B,N}U_n + 2I_N, \quad U_0 = 2I_N.$$

Indeed, the distance between $U_{n+1}$ and $U_n$ controls the error in the operator norm

$$\|A_{B,N}^{*-1} - U_n\| = 2\left\| \sum_{k=n+1}^{+\infty} R_{B,N}^k \right\| \le 2\|R_{B,N}^{n+1}\| \frac{1}{1 - \|R_{B,N}\|} = \frac{\|U_{n+1} - U_n\|}{1 - \|R_{B,N}\|}.$$

Concerning the computation of $\partial_t \mu_S$, we have noticed in Section 3.3, see (3.10), that we should invert

$$\mathcal{A}_{D,N} := A_{S,N}^* - C_{D,N} A_{B,N}^{*-1} D_{D,N},$$

where $C_{D,N}$ and $D_{D,N}$ account for the interactions between the bottom and the free surface. If it has been established that $A_{B,N}^*$ is a perturbation of $\frac{1}{2}I_{N_S}$, it is not the case of $\mathcal{A}_{D,N}$ because of the asymptotic behavior of $C_{D,N}$ and $D_{D,N}$ when the free surface is far from the bottom. Indeed, studying the behavior when $z_S - z_B = iX + o(X)$ for large $X \in \mathbb{R}_+$, we get from the definition of $C_{D,N}$ and $D_{D,N}$ that

$$\begin{aligned} (C_{D,N})_{i,j} &= A_{tw} \frac{L_B}{N_B} \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_B(j)}{L/\pi} \right) z_{B,e}(j) \right] \\ &= A_{tw} \frac{L_B}{N_B} \operatorname{Re} \left[ \frac{1}{2Li} (-iz_{B,e}(j)) \right] + o(1) \\ &= -A_{tw} de_B \frac{\operatorname{Re} z_{B,e}(j)}{2L} + o(1). \end{aligned}$$

and

$$(D_{D,N})_{i,j} = \frac{L_S}{N_S} \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) \right] = de_S \frac{\operatorname{Re} z_{S,e}(j)}{2L} + o(1).$$

Using the decomposition of $A_{S,N}^* = \frac{1}{2}(I_{N_S} - R_{S,N})$ and $A_{B,N}^{*-1} = 2(I_{N_B} + \tilde{R}_{B,N})$ we conclude that

$$\begin{aligned} \mathcal{A}_{D,N} &= \frac{1}{2} I_{N_S} + \frac{A_{tw} de_B de_S}{2L^2} \left( \operatorname{Re} z_{B,e}(j) \right)_{i,j} \left( \operatorname{Re} z_{S,e}(j) \right)_{i,j} + \hat{R}_{B,N} + o(1) \\ &= \frac{1}{2} I_{N_S} + \frac{A_{tw} de_B de_S}{2L^2} \left( \sum_{k=0}^{N_B} \operatorname{Re} z_{B,e}(k) \operatorname{Re} z_{S,e}(j) \right)_{i,j} + \hat{R}_{B,N} + o(1) \\ &= \frac{1}{2} I_{N_S} + \frac{A_{tw} de_S}{2L^2} \left( \operatorname{Re} z_{S,e}(j) \operatorname{Re} \int_0^{L_B} z_{B,e}(e) \, de \right)_{i,j} + \hat{R}_{B,N} + o(1) \end{aligned}$$
$$= \frac{1}{2} \mathrm{I}_{N_S} + \frac{A_{tw} de_S}{2L} \left( \mathrm{Re} \, z_{S,e}(j) \right)_{i,j} + \hat{R}_{B,N} + o(1) = \frac{1}{2} \tilde{\mathcal{A}}_N + \hat{R}_{B,N} + o(1),$$

where $de_B = L_B/N_B$ and $de_S = L_S/N_S$.

We are then first interested in inverting such matrices

$$\tilde{\mathcal{A}} = \mathrm{I}_{N_S} + (a_j)_{i,j}.$$

**Lemma 4.1.** If $1 + \sum_k a_k \neq 0$, then $\tilde{\mathcal{A}}$ is invertible and we have

$$\tilde{\mathcal{A}}^{-1} = \mathrm{I}_{N_S} - \frac{1}{1 + \sum a_k} (a_j)_{i,j}.$$

*Proof.* It is enough to check that the right hand side matrix is the right inverse to $\tilde{\mathcal{A}}$, namely

$$\tilde{\mathcal{A}} \left( \mathrm{I}_{N_S} - \frac{1}{1 + \sum a_k} (a_j)_{i,j} \right) = \mathrm{I}_{N_S} + (a_j)_{i,j} - \frac{1}{1 + \sum a_k} (a_j)_{i,j} - \frac{1}{1 + \sum a_k} (a_j)^2_{i,j} \\ = \mathrm{I}_{N_S} + \left( \frac{a_j(1 + \sum a_k)}{1 + \sum a_k} - \frac{a_j}{1 + \sum a_k} - \frac{a_j \sum a_k}{1 + \sum a_k} \right)_{i,j} = \mathrm{I}_{N_S}.$$

□

In our case, we have

$$1 + \sum a_j = 1 + \frac{A_{tw} de_S}{L} \sum_j \mathrm{Re} \, z_{S,e}(j) \\ = 1 + \frac{A_{tw}}{L} \mathrm{Re} \int_0^{L_S} z_{S,e}(e) \, de + o(1) = 1 + A_{tw} + o(1) \quad (4.1)$$

which is non zero for the single-fluid formulation where $A_{tw} = 1$, but also for the bi-fluid water-waves equations where $A_{tw} > -1$ if $\rho_F > 0$.

By the previous lemma and (4.1), we naturally set

$$\tilde{\mathcal{A}} := \mathrm{I}_{N_S} + \frac{A_{tw} de_S}{L} \left( \mathrm{Re} \, z_{S,e}(j) \right)_{i,j}, \quad \tilde{\mathcal{A}}_{-1} := \mathrm{I}_{N_S} - \frac{A_{tw} de_S}{L(1+A_{tw})} \left( \mathrm{Re} \, z_{S,e}(j) \right)_{i,j}$$

and

$$\mathcal{R} := \mathrm{I}_{N_S} - 2\tilde{\mathcal{A}}_{-1}\mathcal{A}_{D,N},$$

so that

$$\mathcal{A}_{D,N} = \frac{1}{2} \tilde{\mathcal{A}} (\mathrm{I}_{N_S} - \mathcal{R}),$$

where we have neglected the error made in (4.1).

We then have for $\|\mathcal{R}\| < 1$ (because $\mathcal{A}_{D,N} = \frac{1}{2}\tilde{\mathcal{A}} + o(1)$)

$$\partial_t \mu_S = \mathcal{A}_{D,N}^{-1} \text{RHS} = 2(\mathrm{I}_{N_B} - \mathcal{R})^{-1} \tilde{\mathcal{A}}_{-1} \text{RHS} = 2 \sum_{k=0}^{+\infty} \mathcal{R}^k \tilde{\mathcal{A}}_{-1} \text{RHS}.$$
which can be written as a fixed point procedure

$$ 
\begin{aligned}
u_{n+1} &= 2 \sum_{k=0}^{n+1} \mathcal{R}^k \tilde{\mathcal{A}}_{-1} \text{RHS} = \mathcal{R} \left( 2 \sum_{k=0}^{n} \mathcal{R}^k \tilde{\mathcal{A}}_{-1} \text{RHS} \right) + 2 \tilde{\mathcal{A}}_{-1} \text{RHS} \\
&= \mathcal{R} u_n + u_0, \quad u_0 = 2 \tilde{\mathcal{A}}_{-1} \text{RHS}.
\end{aligned}
$$

Concerning the vortex method, we need to inverte $B_{B,N}$ which appears in the determination of $\gamma_B$ in the vortex formulation, see Sections 3.1.1, 3.1.2, 3.1.3. Such an operator is related to the classical vortex method in domains with boundaries, see the operator $B$ in [4, Section 3]. The discrete version

$$ 
\begin{aligned}
B_{B,N}(i,j) &= \frac{L_B}{N_B} \operatorname{Im} \left[ \frac{z_{B,e}(\tilde{e}_{B,i})}{2Li} \cot \left( \frac{z_B(\tilde{e}_{B,i}) - z_B(e_{B,j})}{L/\pi} \right) \right] \quad \forall (i,j) \in [1, N_B - 1] \times [1, N_B], \\
B_{B,N}(N_B, j) &= \frac{L_B}{N_B} \quad \forall j \in [1, N_B],
\end{aligned}
$$

is invertible if $\tilde{e}_{B,i} \in (e_{B,i}, e_{B,i+1})$ are close to a uniform distribution $\tilde{\theta}_i = (\theta_i + \theta_{i+1})/2$. Moreover [4, Theorem 2.1] states that $\gamma_{B,N} = B_{B,N}^{-1} \text{RHS}_{VB,N}$ is a good approximation of $\gamma_B$. Even though it is invertible, the matrix is not well-conditioned and cannot be seen as a Neumann series, except relating $B^{-1}$ to $A^{-1}$ through [4, Equations (3.19) & (3.22)]. However, this step only needs to be performed once at the beginning of the numerical integration (since this matrix does not evolve with time). It is thus worth inverting it accurately.

The last matrix to inverte is

$$ A_{V,N} := A_{S,N} - C_{V,N} B_{B,N}^{-1} D_{V,N} $$

to compute $\partial_t \gamma_S$, see (3.13), which can be also inverted by Neumann series as we did for $A_{D,N}$.

## 5. NUMERICAL RESULTS AND CONVERGENCE

In order to compare and validate the various numerical approaches, we have considered three different initial conditions and test cases when the bottom is flat $\Gamma_B = \mathbb{T}_L \times \{-h_0\}$. In all test cases, we used $N_S = N_B = N$.

While simulations were performed using different values of the Atwood number, we report here simulations with $A_{tw} = 1$, i.e. the single-fluid formulation, which can be compared with existing solutions. We also neglect surface tension in all the test cases below. This makes the tests below more demanding, since the surface tension generally has regularising effects.

The parameter $\alpha$ was set to $\alpha = 1$, implying that the points are advected tangentially at the velocity of the lower fluid. Also all test cases presented here include a flat bottom (though the cases of infinite depth or variable bottom can also be handle using the same code). The simulations are made non-dimensional, with a length-scale based on the water depth, such that $h_0 = 1$, and a time-scale based on gravity, such that $g = 1$.
We will numerically investigate the stability of our scheme in a few test cases. We will then validate the numerically obtained solutions against analytical solutions, when available. When not available, we will simply check convergence to an N independent solution in the limit of large N. Another useful validation concerns conserved quantities. The first thing to be checked is the conservation of mass

$$ M(t) := \rho_F \mathrm{Vol} \mathcal{D}_F(t) = \rho_F \iint_{\mathcal{D}_F(t)} \mathrm{div} \begin{pmatrix} 0 \\ x_2 \end{pmatrix} \mathrm{d}\mathbf{x} = \rho_F \int_{\partial \mathcal{D}_F(t)} \begin{pmatrix} 0 \\ x_2 \end{pmatrix} \cdot \tilde{\mathbf{n}} \mathrm{d}\sigma(\mathbf{x}) \\ = \rho_F \int_0^{L_S} z_{S,2}(e) z_{S,e,1}(e) \mathrm{d}e - \rho_F \int_0^{L_B} z_{B,2}(e) z_{B,e,1}(e) \mathrm{d}e $$

and the total energy

$$ E(t) := \frac{1}{2} \iint_{\mathcal{D}_F(t)} \rho_F |\mathbf{u}|^2 + \iint_{\mathcal{D}_F(t)} \rho_F g x_2 \\ = \frac{\rho_F}{2} \iint_{\mathcal{D}_F(t)} \nabla \phi_F \cdot \nabla^\perp \psi + \frac{\rho_F g}{2} \iint_{\mathcal{D}_F(t)} \mathrm{div} \begin{pmatrix} 0 \\ x_2^2 \end{pmatrix} \mathrm{d}\mathbf{x} \\ = -\frac{\rho_F}{2} \int_0^{L_S} \mathrm{Re} \left[ \hat{\mathbf{u}}_F(z_S(e)) z_{S,e}(e) \right] \psi(z_S(e)) \mathrm{d}e \\ + \frac{\rho_F}{2} \int_0^{L_B} \mathrm{Re} \left[ \hat{\mathbf{u}}_F(z_B(e)) z_{B,e}(e) \right] \psi(z_B(e)) \mathrm{d}e \\ + \frac{\rho_F g}{2} \int_0^{L_S} z_{S,2}^2(e) z_{S,e,1}(e) \mathrm{d}e - \frac{\rho_F g}{2} \int_0^{L_B} z_{B,2}^2(e) z_{B,e,1}(e) \mathrm{d}e $$

where the expression of $\hat{\mathbf{u}}_F$ and $\psi$, given in (2.14) and (2.22), has to be considered with the limit formulae of Appendix A.

We also successfully checked conserved integrals for domains with horizontal symmetries [1] (not detailed here).

5.1. **Case 1: linear water waves.** The first test case we consider is that of simple waves of small amplitude, it takes the form provided by Stokes at first order

$$ \eta(t, x) = A \cos(kx - \omega t), \quad \Phi(t, x, y) = A \frac{\omega}{k} \frac{\cosh k(y + h_0)}{\sinh kh_0} \sin(kx - \omega t), \quad (5.1) \\ \text{with } \omega = \sqrt{gk \tanh kh_0}, \text{ and } u_x = \partial_x \Phi, u_y = \partial_y \Phi. $$

Our initial condition is thus

$$ \eta_0(x) = A \cos(kx), \quad \mathbf{u} \cdot \mathbf{n} = A \sin kx \sqrt{gk \tanh kh_0}. \quad (5.2) $$

We also consider Stokes waves at second order, where the initial data is slightly modified (see [24]):

$$ \eta_0(x) = A \cos(kx) + kA^2 \frac{3 - \tanh^2 kh_0}{4 \tanh^3 kh_0} \cos(2kx), \\ \mathbf{u} \cdot \mathbf{n} = \frac{A \sin kx \sqrt{gk \tanh kh_0}}{\sqrt{1 + k^2 A^2 \sin^2(kx)}} \left( 1 + \frac{kA}{\tanh kh_0} \cos(kx) \right). $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Convergence of the numerical integration after 10 revolution periods, for simple waves with amplitude (a) 0.01 and (b) 0.001. The error (compared to the initial condition) defined as $\mathcal{E} = \|z(t_f) - z(0)\|_{\infty}$ is represented as a function of the number $N$ of points used in the numerics. The black dashed curve corresponds to the first order simple wave, whereas the solid red curve corresponds to the second order Stokes solution. The vortex and the dipole methods are indistinguishable in these graphs.</figcaption></figure>

In our simulations, we used $\mathbf{u} \cdot \mathbf{n}$ as boundary condition to numerically find $\gamma_0$ and $\mu_0$ (see § 3.1.1 and § 3.1.3). This contrasts with [6] which provides the analytical expression for both $\gamma$ and $\mu$ in the limit of a vanishing amplitude. This distinction is small at this stage, but turns out to be important at later stage (see our third test case below).

We consider numerically a wave number $k = 1$, and a domain of horizontal extent $L = 2\pi$. We integrate our simulations for a time $t_f = 10k/\omega$. We introduce a CFL-number

$$CFL = \min(|\dot{z}_S| \Delta t / (|z_{S,e}| \text{de})),$$

based on the Lagrangian velocity of the points sampling the interface. We have checked that our numerical code is stable for a CFL-number $CFL \le CFL^*$, we observed numerically that $CFL^* \in [0.25, 0.5)$. We can then test the convergence of our numerical solution by varying the spatial resolution $N$. In order to keep temporal truncation error, we then used $CFL \le 1/10$. Convergence is then assessed by comparing the final state to the initial condition $\mathcal{E} = \|z(t_f) - z(0)\|_{\infty}$. The resulting errors for waves of amplitudes $A = 10^{-2}$ and $A = 10^{-3}$ respectively are represented in Figure 2(a,b). The error, defined in $L^\infty$-norm, decreases with increasing resolution to a value which decreases with the amplitude of the wave. This can be interpreted as the signature of weak non-linear corrections. Interestingly, the vortex and the dipole methods are indistinguishable in these graphs.

5.2. **Case 2: solitary waves.** Our second test case concerns solitary waves. We consider here an extension of solitary waves to a periodic domain. This involves the cnoidal wave solution for the Green-Naghdi equation.
To present this solution, first we define the Jacobi elliptic functions $cn(u, m) := \cos \varphi(u, m)$ as the inverse of the elliptic integral

$$u(\varphi, m) := \int_{0}^{\varphi} \frac{d\theta}{\sqrt{1 - m \sin^2 \theta}}$$

where $m \in (0, 1)$. We also need the complete elliptic integral of the first and second kind:

$$K(m) := \int_{0}^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1 - m \sin^2 \theta}} \quad \text{and} \quad E(m) := \int_{0}^{\frac{\pi}{2}} \sqrt{1 - m \sin^2 \theta} \, d\theta.$$

For a given period $L$, amplitude $A$ and depth $h_0$, we determine the nonlinearity parameter $m \in (0, 1)$ verifying the following dispersion relation

$$AL^2 = \frac{16}{3} m K^2(m) \frac{h_0^2}{g} c^2(m),$$

where the velocity $c$ is given by

$$c(m) := \sqrt{gh_0\left(1 + \frac{\eta_1(m)}{h_0}\right)\left(1 + \frac{\eta_2(m)}{h_0}\right)\left(1 + \frac{\eta_3(m)}{h_0}\right)}$$

with

$$\eta_1(m) := -\frac{A}{m} \frac{E(m)}{K(m)}, \quad \eta_2(m) := \frac{A}{m} \left(1 - m - \frac{E(m)}{K(m)}\right), \quad \eta_3(m) := \frac{A}{m} \left(1 - \frac{E(m)}{K(m)}\right).$$

Once $m$ is obtained, the Green-Naghdi soliton is then defined by the surface elevation

$$\eta(t, x) = \eta_2 + A \operatorname{cn}^2 \left( \frac{2K(m)}{L} (x - ct), m \right).$$

As it translates to the right with velocity $c$, we write

$$\mathbf{u} \cdot \mathbf{n}|_{(0,x,\eta(0,x))} = -c \frac{\partial_x \eta(0, x)}{\sqrt{1 + |\partial_x \eta(0, x)|^2}},$$

which uniquely defines $\gamma_0$ and $\mu_0$, see § 3.1.1 and § 3.1.3.

For more details on cnoidal waves, we refer to [12, 23], and references therein.

We consider numerically a domain of large extent to allow for a localized solitary wave, and choose $L = 40\pi$. For an amplitude $A = 0.1$, we computed a circulation $\gamma = 3 \cdot 10^{-3}$. We checked numerically that setting $\gamma = 0$ did not alter the solution in this case (i.e. the circulation is weak enough not to affect the solution). The results are presented in Figure 3. Figure 3(a) highlights the modification induced on the above initial condition after one full period, i.e. $t^* = L/c \simeq 124.23$. The difference (represented in dashed blue) appears dominated by a correction on the velocity $c$, presumably due to higher order correction terms. We performed simulations up to $t = 10t^*$ which confirmed that no change of shape is observed, but the phase lag with the analytical solution increases with time.

Both methods appear stable in time and no significant change on the numerical solution was observed when the time-step was reduced by a factor of 2.

The numerical convergence of both methods is demonstrated in Figure 3(b). Here the solitons are localized and a small difference of phase velocity yields a large $L^\infty$-norm, not
<figure><img src="image_5.png" /><figcaption>FIGURE 3. (a) The initial solitary wave (dashed black) and the numerical integration with $N = 2048$ using the dipole formulation up to $t = 124.23$, i.e. after one revolution of the analytical solution. (b) Convergence of the dipole method (blue) and vortex method (red). The error, measured using the Hausdorff distance, for a given resolution $N$ is computed relative to $N = 2048$ for the same method at time $t = 124.23$. The black dotted line illustrates second order convergence. The full computational domain extends to $L = 40\pi$.</figcaption></figure>

reflecting the actual distance between the two curves. We thus use a discrete Hausdorff distance between two solutions to measure the error (see Fig. 3.b):

$$ \mathcal{E} = d_H(z_1, z_2) \equiv \max \left( \max_i(\min_j(|z_1(i) - z_2(j)|)), \max_j(\min_i(|z_1(i) - z_2(j)|)) \right), \quad (5.3) $$

where both curves were interpolated using $2^{17} \simeq 130.000$ points in order to identify accurately enough the closer points on both curves.

Both methods exhibit an approximately second order convergence in space (as highlighted by the black dotted line). Both the vortex and the dipole methods yield very good results on this test case as well.

For this rather long integration (up to $t = t^* \simeq 124.23$), we observe that the volume is conserved up to fluctuations of the order of $10^{-8}$ and the total energy (kinetic and potential) is very weakly dissipated (or the order of $10^{-4}$ over the full period).

### 5.3. Case 3: wave breaking.
Finally we turn to a strongly non-linear problem, that of a wave breaking. We consider a flat bottom and a mean water depth $h = 1$. We consider an initial interface of the form

$$ \eta_0(x) = A \cos(kx) \quad \text{with } A = \frac{1}{2} \text{ and } k = 1. \quad (5.4) $$

Our initial velocity stems from (5.2) but is now used with large amplitude. In doing so it is important to relax the approximation $n \simeq e_y$. We start with (5.1) which yields the
<table><caption>TABLE 1. Final integration time before the occurrence of a numerical instability for various numerical schemes (the F-vortex method is introduced in § 6.1). We also integrated the OEC dipole formulation (introduced later in the paper) with 4096 points and obtained a stable solution (nearly undistinguishable from the 2048 points simulation) until $t = 3.61$.</caption><thead><tr><th></th><th>Vortex</th><th>Dipole</th><th>OEC-dipole</th><th>F-vortex</th></tr></thead><tbody><tr><td>256</td><td>2.40</td><td>3.02</td><td>3.06</td><td>2.81</td></tr><tr><td>512</td><td>2.03</td><td>3.04</td><td>3.37</td><td>2.80</td></tr><tr><td>1024</td><td>1.73</td><td>3.13</td><td>3.39</td><td>2.86</td></tr><tr><td>2048</td><td>1.15</td><td>3.10</td><td>3.60</td><td>2.93</td></tr></tbody></table>

initial condition on velocity

$$u_x = \partial_x \Phi = A \sqrt{\frac{gk}{\tanh kh}} \cos kx, \quad u_y = \partial_y \Phi = A \sqrt{gk \tanh kh} \sin kx.$$

Using expression of the normal

$$\mathbf{n} = \frac{1}{\sqrt{1 + (\partial_x \eta)^2}} \begin{pmatrix} -\partial_x \eta \\ 1 \end{pmatrix} = \frac{1}{\sqrt{1 + k^2 A^2 \sin^2 kx}} \begin{pmatrix} kA \sin kx \\ 1 \end{pmatrix},$$

we get

$$\mathbf{u} \cdot \mathbf{n} = \frac{A \sin kx \sqrt{gk \tanh kh}}{\sqrt{1 + k^2 A^2 \sin^2 kx}} \left( 1 + kA \frac{1}{\tanh kh} \cos kx \right) \quad (5.5)$$

which resumes to (5.2) in the limit of vanishing amplitudes. It is interesting to note that this differs from [6], which used the $\gamma_0$ and $\mu_0$ constructed from the small amplitude approximation.

It should also be stressed that the construction of a sensible initial condition at large amplitude (wave breaking) from the simple wave analytics is necessarily based on simplifying assumption as there is no analytical solution in this strongly non-linear limit.

This situation is in a fully non-linear regime and more challenging numerically than the two previous test cases.

We observed here, as reported by [6], that for fully-nonlinear configurations the vortex method is unstable to a high wave number instability. As the resolution is increased and higher wave numbers are resolved, the integration time before an instability occurs decreases (see first column in Table 1).

Such is not the case (again as stressed by [6]) for the dipole method, for which the integration time is fairly independent of the resolution (see second column in Table 1).

As we have shown in the first two cases, both the vortex method and the dipole method are stable and can be used in situations involving a small to moderate curvature. In the case of huge curvature (such as wave-breaking), the vortex method becomes impractical and does not converge (as the stability decreases with increasing resolution).
For both methods, we observed that the volume and total energy are conserved up to less than 1%.

The case of the dipole method is not quite as severe, but we observe a numerical instability at a time independent of $N$, well before the splash occurs. The instability of the dipole method appears to develop first in the form of points approaching each other in the direction tangent to the interface, in a manner similar to the so called ‘phantom traffic jam’ instability. In order to delay the formation of this instability, we introduced an ‘odd-even coupling’ (OEC) at the end of each time-step in the form of an arbitrary regularization on $\mu$ by replacing at the end of each time step, $\partial_t\mu_S(e_k)$ with $(\partial_t\mu_S(e_{k-1}) + 2\partial_t\mu_S(e_k) + \partial_t\mu_S(e_{k+1})) / 4$. Using this approach we could extend the integration time by nearly 10% (see Table 1).

For large grids, this coupling procedure does not alter the simulations on the first two test-cases such as simple waves or solitary waves, but does stabilise the numerical integration of wave breaking.

The OEC approach introduces a stabilisation, which however does not affect the overall convergence of the scheme. At times $t = 1$, $t = 2$ and $t = 3$, the convergence of the OEC approach is undistinguishable from that of the dipole method, see figure 4. Besides the OEC vanishes continuously in the limit of large grids.

It is worth stressing that contrarily to other regularization techniques previously used on this problem (and discussed in the following section), the above OEC does not involve any arbitrary small parameter $\varepsilon$ other than the grid space. This approach is thus free of the risk to present results with vanishing distance between points and yet a finite regularization.

The curves marking the interface between water and air, are generally not graphs for this test-case. We thus stick to the discrete Hausdorff distance between two curves, see (5.3), to measure the error (see Figs. 4 and 5).

It is worth stressing that the very same test case has been recently investigated using the Navier-Stokes equations and a Finite Element discretization, and that convergence has been achieved to the solution portrayed on Fig. 5.a as the Reynolds number is increased [32].

## 6. COMPARISON WITH REGULARIZATION STRATEGIES

### 6.1. Fourier filtering for the vortex method.
Since the Vortex methods appears unstable (see [6] and Table 1), some Fourier filtering can be introduced. The strategy of a filtered Vortex method is for example followed by [3].

The filtering corresponds to a product in Fourier space of both $z_S$ and $\gamma_S$ with a filter function. We used here the filter function introduced in [7] (see eq. (2.14) in this reference)

$$ \hat{F} = \frac{1}{2} - \frac{1}{2} \tanh \left( \frac{2|k|\pi/N_S - \xi_0}{d} \right). $$

Two parameters are introduced in the above. $\xi_0$ locates the center of the transition zone (usually some fraction of $\pi$) and $d$ controls the width of the transition zone. Following [7], we used $\xi_0 = \pi/4$ and $d = \pi/40$.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Convergence of the dipole method (a) and of the OEC-dipole (b). In both cases the reference solution is that obtained with the dipole method for $N = 2048$. The error, measured by the Hausdorff distance, is represented at time $t = 1.00$ (solid blue), $t = 2.00$ (dashed red) and time $t = 3.00$ (dot dashed green). The dotted black curve indicates second order decrease.</figcaption></figure>

The Fourier filtered method (referred to as ‘F-vortex’ in the text) yields an increased stability and thus longer integration. Remarkably the final time of integration for the F-vortex method is independent of $N$ (see Table 1) and the method does converge see Fig. 6(a). Although the simulation has been extended in time far beyond the unfiltered vortex method, the solution does converge (to first order) to that of the dipole method, see Fig. 6(b). The final integration time however remains shorter than for the dipole method (let alone the OEC-dipole method).

The use of $\xi_0 = \pi/4$ in the above tests (guided by [7]) yields an ‘effective’ resolution of approximately $N/4$, (though with a larger stability than the pure vortex method with $N/4$). Increasing the cut-off frequency, say to $\xi_0 = \pi/2$ instead of $\xi_0 = \pi/4$ yields a less stable scheme. The observed time for instability with $\xi_0 = \pi/2$ was $t = 2.64$ for $N = 512$ and $t = 2.24$ for $N = 1024$.

We should also stress that [7] introduced, in the case of a very stiff initial data, a filtering on the dipole method. This interesting approach stabilizes the dipole method, thus allowing for longer time integration.

## 6.2. Curve-offset method.
Another approach to regularise the boundary integral consists in considering that the vortices are located at a finite distance above the free surface (e.g. [13, 34, 35, 33]). This finite offset prevents any singularity in the integration as the vortices are fictiously located at

$$ (X_j, Y_j) = (x_j, y_j) + L_d \mathbf{n} \quad \text{with } L_d = \delta L/N. $$

While no singularity can occur with this technique, the kernel in the vortex formula- tion takes the form cot($\pi(x_{s,i} - x_{v,j})/L$) where $(x_{s,i})$ are located on the free surface
<figure><img src="image_7.png" /><figcaption>FIGURE 5. (a) Time evolutions of the breaking wave with $N = 4096$ and the OEC-dipole discretization from the $\frac{1}{2} \cos(x)$ profile at $t = 0$ (see (5.4)-(5.5)) to the splash at $t \approx 3.6$. (b) Plots at time $t = 1, t = 2, t = 3, t = 3.6$. (c) Solutions using the OEC-dipole discretization with 256 points (dashed blue) and 4096 points (solid red) at time $t=3.0$. (d) Convergence of the OEC-dipole toward the OEC-dipole $N = 4096$ curve used as reference at time $t = 1.0$ (solid blue), $t = 2.0$ (dashed red) and time $t = 3.0$ (dot dashed green). The dotted line indicates the exact second order convergence slope.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Convergence of the F-vortex method at time $t = 2.75$; (a) compared to $N = 2048$ with the same method (used as a reference). (b) compared to the dipole method with $N = 2048$. The dotted black line indicates second order convergence.</figcaption></figure>

whereas $(x_{v,i})$ are located a finite distance above the free surface. It is finite but of order $\cot(\pi L_d / L)$, which becomes large when $L_d$ is small.

This method is largely discussed in [33, p. 928] where the author introduces $L_d$, but also a second regularization which consists in replacing the Green kernel $G(x_{s,i}, x_{v,j})$, behaving as $\ln|x_{s,i} - x_{v,j}|$ (see (2.1)), with $\ln|x_{s,i} - x_{v,j}| + b$ where $b$ is not necessarily small ($b \in [0; 10\,000]$). The above $b$ term is difficult to justify from a mathematical point of view, and not quite a small perturbation.

Instead of the above procedure, we consider a regularization inspired from the vortex-blob method. We thus introduce a second regularizing parameter $\epsilon_N > 0$ replacing the previous kernel functions by $\cot(\pi(x_{s,i} - x_{v,j})/L + \epsilon_N n)$.

In practice, when testing this approach, we have considered the $L_d$ regularization (finite distance of the vortices from the interface) and a parameter $\epsilon_N$ on the form of the blob method (i.e. regularizing the kernel). Following [33], we explicitely relate $L_d$ to $1/N$. Also $\epsilon_N$ is taken to vanish as $1/N$ to try to assess the convergence properties of this approach.

We present in Fig. 7 the numerical simulations performed with the curve-offset method applied to the vortex formulation. We used (5.2) as initial condition. This configuration corresponds to the wave breaking and the results should be compared with those of Fig. 5.

We considered four different regularizations weights, namely $L_d = L/N \epsilon_N = 1/2N$, $L_d = 2L/N \epsilon_N = 1/2N$, $L_d = L/2N \epsilon_N = 1/2N$, $L_d = L/N \epsilon_N = 1/4N$. We report the various numerical solutions at $t = 3.68$ in Fig. 7(a,b). We first note that all curves are significantly different from the results presented in Fig. 5. The wave is much slower with the regularization. The unregularized method, for which we observed convergence toward the Euler solution on various test cases, has already reached splashing at that time. The shape of the obtained numerical solution also strongly depends on $L_d$, but only weakly on $\epsilon_N$ (the green and black curves being extremely close to another). We performed further tests, which confirmed the weak influence of the $\epsilon_N$ regularization on the numerical solution.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Simulations of wave breaking with initial condition (5.2) using the curve-offset method. In graphs (a) and (b), profiles of the interface obtained with $N = 512$ at time $t = 3.68$ with $L_d = L/N$, $\epsilon_N = 1/2N$ (black), with $L_d = 2L/N$, $\epsilon_N = 1/2N$ (blue), with $L_d = L/2N$, $\epsilon_N = 1/2N$ (red), and with $L_d = L/N$, $\epsilon_N = 1/4N$ (green). (c) Evolution of the Hausdorff error, between $N=2048$ and $N=256, 512, 1024$, at time $t = 1.0$ (solid blue), $t = 2.0$ (dashed red) and time $t = 3.0$ (dot dashed green). The dotted line indicates an ideal $N^{-2}$ scaling.</figcaption></figure>

A puzzling property of this approach is that, for a given choice of regularization, say $L_d = L/N$, $\epsilon_N = 1/2N$, some convergence is achieved as $N$ is increased (see Fig. 7c). Yet the numerical curve then converges toward a curve which could seem plausible, but significantly differs from the unregularized solution. A final concern with this approach is that the total energy of the wave is not conserved (and varies by $\sim 30\%$ through the simulation).

## 7. DISCUSSION

We derived a numerical strategy to discretize inviscid water waves in the case of over-
turning interfaces (i.e. when the water-air interface is not a graph). We showed that this
discretisation can be used up to the splash (i.e. when the interface self intersects). No
filtering or regularisation were introduced other than numerical discretisation. In the most
severe case of a splashing wave, an odd-even coupling was introduced. It vanishes in the limit of a large number of points.

This formulation opens the way for further studies. In particular, we want to study the possibility of a finite-time singularity formation at the tip of a breaking wave. No singularity was observed with the initial condition considered here. We also want to investigate the effect of an abrupt jump in water height. Finally, the triple interface of water, air and land (i.e. the sloping beach problem) still needs to be addressed.

## ACKNOWLEDGEMENTS

The authors are grateful to David Lannes and Alan Riquier for discussions. This work was partially supported by the ANR project 'SINGFLOWS' (ANR-18-CE40-0027-01), the IMPT project 'Ocean waves', and the CNES-Tosca project 'Maeva'.

## APPENDIX A. COTANGENT KERNEL

We begin this appendix by relating the cotangent kernel and the usual kernel in $\mathbb{R}^2$. Let $f$ a $L$-periodic real function and $z$ a curve verifying $z(e + L) = z(e) + L$, then (we recall the notation $\hat{u} = u_1 - iu_2$,

$$ 
\begin{aligned}
\widehat{K}_{\mathbb{R}^2}[|z_e|^{-1}f](x) &= \frac{1}{2\pi i} \int_{-\infty}^{\infty} \frac{1}{x-z(e')} f(e') \, de' = \frac{1}{2\pi i} \int_{0}^{L} \sum_{k=-\infty}^{+\infty} \frac{1}{x-z(e') - Lk} f(e') \, de' \\
&= \frac{1}{2\pi Li} \int_{0}^{L} \sum_{k=-\infty}^{+\infty} \frac{1}{\frac{1}{L}(x-z(e')) - k} f(e') \, de' \\
&= \frac{1}{2Li} \int_{0}^{L} \cot\left(\frac{x-z(e')}{L/\pi}\right) f(e') \, de'
\end{aligned}
$$

because it is well known for $x \in (0,1)$ that

$$ \pi \cot(\pi x) = \lim_{N \to \infty} \sum_{k=-N}^{N} \frac{1}{x+n} = \frac{1}{x} + \lim_{N \to \infty} \sum_{k=1}^{N} \frac{2x}{x^2-n^2} $$

which can be extended to $C \setminus \mathbb{Z}$ by unicity of the analytic extension. This computation can be found in several formal derivations to get the Biot-Savart formula without using the Green kernel in $T_L \times \mathbb{R}$, see for instance [3].

The limit of Cauchy integrals from above provides

$$ 
\begin{aligned}
\lim_{s \to 0^+} \frac{1}{2Li} \int_0^L \cot \left( \frac{x - z(e')}{L/\pi} \right) f(e') \, de' \Big|_{x=z(e)+sn} \\
= \frac{1}{2Li} \operatorname{pv} \int_0^L \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) f(e') \, de' - \frac{1}{2} \frac{f(e)}{z_e(e)}
\end{aligned}
$$

and from below we have

$$ \lim_{s \to 0^-} \frac{1}{2Li} \int_0^L \cot \left( \frac{x - z(e')}{L/\pi} \right) f(e') \, de' \Big|_{x=z(e)+sn} $$
$$ = \frac{1}{2Li} \mathrm{pv} \int_{0}^{L} \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) f(e') \, de' + \frac{1}{2} \frac{f(e)}{z_{e}(e)}. $$

From this formula, we recover that the tangential component $\mathbf{u} \cdot \tau = \mathrm{Re}(\hat{\mathbf{u}}\tau)$ has a jump, whereas the normal component $\mathbf{u} \cdot \mathbf{n} = -\mathrm{Im}(\hat{\mathbf{u}}\tau)$ is continuous.

We can also note that $\int_0^L \mathrm{Re} \left[ \frac{z_e(e)}{2Li} \cot \left( \frac{z(e)-z(e')}{L/\pi} \right) \right] f(e') \, de'$ is actually a classical integral whereas $\mathrm{pv} \int_0^L \mathrm{Im} \left[ \frac{z_e(e)}{2Li} \cot \left( \frac{z(e)-z(e')}{L/\pi} \right) \right] f(e') \, de'$ only makes sense in terms of principal value.

As usual concerning desingularization of the principal value, we first note that

$$ \begin{aligned} \mathrm{pv} \int_0^L z_e(e') \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) \, de' &= \frac{L}{\pi} \mathrm{pv} \int_{-\infty}^{\infty} \frac{z_e(e')}{z(e) - z(e')} \, de' \quad (A.1) \\ &= -\frac{L}{\pi} \lim_{\epsilon \to 0^+} \left( \left[ \ln(z(e) - z(e')) \right]_{-\infty}^{e-\epsilon} + \left[ \ln(z(e) - z(e')) \right]_{e+\epsilon}^{+\infty} \right) \\ &= -\frac{L}{\pi} \lim_{\epsilon \to 0^+} \lim_{S \to \infty} \left( \left( \ln(\epsilon\rho) + i\theta \right) - \ln S + \left( \ln S + i\pi - \ln(\epsilon\rho) + i\theta + i\pi \right) \right) = 0 \end{aligned} $$

where $z(e) - z(e \pm \epsilon) = - \pm \epsilon z_e(e) + \mathcal{O}(\epsilon^2) = - \pm \epsilon \rho e^{i\theta} + \mathcal{O}(\epsilon^2)$, hence we can always write

$$ \mathrm{pv} \int \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) f(e') \, de' = \int \cot \left( \frac{z(e) - z(e')}{L/\pi} \right) \frac{f(e')z_e(e) - f(e)z_e(e')}{z_e(e)} \, de' $$

which is now a classical integral of a continuous function.

For more details on singular integrals, we refer to [29, 25, 17], see [4, Sect. 3] for a brief summary.

## APPENDIX B. DISCRETE OPERATORS FOR THE DIPOLE FORMULATION

For $z_B(i) := z_B(e_{B,i})$, $z_S(i) := z_S(e_{S,i})$ and $\mu_S(i) := \mu_S(e_{S,i})$ given, we construct the matrix $A_{B,N}^*$

$$ A_{B,N}^*(i,j) = \frac{L_B}{N_B} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_B(j)}{L/\pi} \right) z_{B,e}(j) \right] \quad \forall i \neq j \in [1, N_B] \times [1, N_B], $$

$$ A_{B,N}^*(i,i) = \frac{1}{2} - \frac{L_B}{N_B} \mathrm{Re} \left[ \frac{1}{4\pi i} \frac{z_{B,ee}(i)}{z_{B,e}(i)} \right] \quad \forall i \in [1, N_S], $$

and $F_{D,N}$

$$ F_{D,N}(i) = - \sum_{j=1}^{N_S} \frac{L_S}{N_S} \mu_S(j) \mathrm{Re} \left[ \frac{z_{S,e}(j)}{2Li} \cot \left( \frac{z_B(i) - z_S(j)}{L/\pi} \right) \right] \quad \forall i \in [1, N_B]. $$

We set $\mu_B = (A_{B,N}^*)^{-1} F_{D,N}$. This operation corresponds to (3.2).
Next, we compute $\gamma_B = \partial_e \mu_B(e)$, $\gamma_S = \partial_e \mu_S(e)$, next $\partial_t z_S$, $\overline{\nabla \phi_F}(z_S(e))$ and $\overline{\nabla \phi_A}(z_S(e))$ where we could extend in the integral on $\Gamma_S$ for $e' = e$ by $\frac{\gamma_S(e)z_{S,ee}(e) - \gamma_{S,e}(e)z_{S,e}(e)}{2\pi i z_{S,e}^2(e)}$.

We compute $A_{S,N}^*$ as

$$A_{S,N}^*(i,j) = A_{tw} \frac{L_S}{N_S} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) \right] \quad \forall i \neq j \in [1, N_S] \times [1, N_S],$$

$$A_{S,N}^*(i,i) = \frac{1}{2} - A_{tw} \frac{L_S}{N_S} \sum_{j \neq i} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) \right] \quad \forall i \in [1, N_S],$$

next

$$C_{D,N}(i,j) = A_{tw} \frac{L_B}{N_B} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_B(j)}{L/\pi} \right) z_{B,e}(j) \right] \quad \forall (i,j) \in [1, N_S] \times [1, N_B].$$

and finally

$$D_{D,N}(i,j) = \frac{L_S}{N_S} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) \right] \quad \forall (i,j) \in [1, N_B] \times [1, N_S].$$

Concerning the right hand side term, we compute

$$G_{D,1,N}(i) = -A_{tw} \sum_{j \neq i} \frac{L_S}{N_S} (\mu_S(i) - \mu_S(j)) \mathrm{Re} \left[ \frac{\pi}{2L^2i} \sin^{-2} \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) (\partial_t z_S(i) - \partial_t z_S(j)) z_{S,e}(j) \right] \\ -A_{tw} \sum_{j \neq i} \frac{L_S}{N_S} (\mu_S(j) - \mu_S(i)) \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) \partial_t z_{S,e}(j) \right] \\ +A_{tw} \sum_{j=1}^{N_B} \frac{L_B}{N_B} \mu_B(j) \mathrm{Re} \left[ \frac{\pi}{2L^2i} \sin^{-2} \left( \frac{z_S(i) - z_B(j)}{L/\pi} \right) \partial_t z_S(i) z_{B,e}(j) \right] \quad (B.1) \\ +\frac{1}{2} \mathrm{Re} \left[ \partial_t z_S(i) \left( (A_{tw}+1) \overline{\nabla \phi_F}(z_S(i)) + (A_{tw}-1) \overline{\nabla \phi_A}(z_S(i)) \right) \right] \\ -\frac{1}{4} \left( (A_{tw}+1) |\overline{(\nabla \phi_F + \mathbf{u}_\gamma)}(z_S(i))|^2 + (A_{tw}-1) |\overline{(\nabla \phi_A(z_S(i)))|^2 \right) \\ -\frac{(A_{tw}+1)\sigma}{2\rho_F} \kappa(z_S(i)) - gA_{tw} \operatorname{Im} z_S(i),$$

for all $i \in [1, N_S]$, whereas

$$G_{D,2,N}(i) = - \sum_{j=1}^{N_S} \frac{L_S}{N_S} \mu_S(j) \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_S(j)}{L/\pi} \right) \partial_t z_{S,e}(j) \right]$$
$$ -\sum_{j=1}^{N_S} \frac{L_S}{N_S} \mu_S(j) \text{Re} \left[ \frac{\pi}{2L^2i} \sin^{-2} \left( \frac{z_B(i) - z_S(j)}{L/\pi} \right) \partial_t z_S(j) z_{S,e}(j) \right] $$

for all $i \in [1, N_B]$.

Remark B.1. It could seem strange that the diagonal terms in $A_{S,N}^*$ are of a different nature than those in $A_{B,N}^*$. It is in fact the same, because we can make use of Appendix A to rewrite (3.2) in the form

$$ \frac{1}{2}\mu_B(e) + \int_0^{L_B} (\mu_B(e') - \mu_B(e)) \text{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - z_B(e')}{L/\pi} \right) z_{B,e}(e') \right] \mathrm{d}e' \\ = - \int_0^{L_S} \mu_S(e') \text{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) z_{S,e}(e') \right] \mathrm{d}e', $$

for which we would defined

$$ A_{B,N}^*(i, i) = \frac{1}{2} - \frac{L_B}{N_B} \sum_{j \neq i} \text{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_B(j)}{L/\pi} \right) z_{B,e}(j) \right] \quad \forall i \in [1, N_S] $$

and we recover the same expression by the discretization of desingularization rule (A.1)

$$ \frac{L_B}{N_B} \sum_{j \neq i} \text{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_B(i) - z_B(j)}{L/\pi} \right) z_{B,e}(j) \right] - \frac{L_B}{N_B} \text{Re} \left[ \frac{1}{4\pi i} \frac{z_{B,e}(i)}{z_{B,e}(i)} \right] = 0. $$

We can therefore use either of these formulations, but it is more interesting to avoid the second derivative $z_{S,ee}$, which tends to destabilize the numerical code when the curvature of the interface becomes large. As the bottom boundary does not depend on time, we can retain the expression in terms of $z_{B,ee}$.

This relation could be also used in the extension for $\partial_t z_S$ mentioned above replacing

$$ \frac{L_S \gamma_S(i) z_{S,ee}(i) - \gamma_{S,e}(i) z_{S,e}(i)}{N_S 2\pi i z_{S,e}^2(i)} \quad \text{with} $$

$$ \frac{L_S \gamma_S(i)}{N_S z_{S,e}(i)} \sum_{j \neq i} \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) - \frac{L_S \gamma_{S,e}(i)}{N_S 2\pi i z_{S,e}(i)}. $$

Unfortunately, this does not improve the stability of the code. The method explained in § 4.2 with the shifted grids in space allows us to avoid $\gamma_{S,e}$, which would appear by the extension by continuity.

## APPENDIX C. DISCRETE OPERATORS FOR THE VORTEX FORMULATION

First, we give the precise equation for the vortex formulation, then we give the discret version of the operators.

We compute

$$ \frac{1}{2} \partial_t \Psi_S(e) = \int_0^{L_S} \text{Re} \left[ \frac{\partial_t \gamma_S(e') z_{S,e}(e) - \partial_t \gamma_S(e) z_{S,e}(e')}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{d}e' $$
$$
\begin{align*}
&+ \int_0^{L_B} \partial_t \gamma_B(e') \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) z_{S,e}(e) \right] \mathrm{d}e' \\
&- \int_0^{L_S} \operatorname{Re} \left[ \pi \frac{\gamma_S(e') z_{S,e}(e) - \gamma_S(e) z_{S,e}(e')}{2L^2 i} \sin^{-2} \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) (\partial_t z_S(e) - \partial_t z_S(e')) \right] \mathrm{d}e' \\
&+ \int_0^{L_S} \operatorname{Re} \left[ \frac{\gamma_S(e') \partial_t z_{S,e}(e) - \gamma_S(e) \partial_t z_{S,e}(e')}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{d}e' \\
&+ \int_0^{L_B} \gamma_B(e') \operatorname{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \partial_t z_{S,e}(e) \right] \mathrm{d}e' \\
&- \int_0^{L_B} \gamma_B(e') \operatorname{Re} \left[ \frac{\pi}{2L^2 i} \sin^{-2} \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \partial_t z_S(e) z_{S,e}(e) \right] \mathrm{d}e' \\
&+ \sum_{j=1}^{N_v} \gamma_{v,j} \operatorname{Re} \left[ \frac{\partial_t z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_{v,j}}{L/\pi} \right) - \frac{\pi z_{S,e}(e) \partial_t z_{S,e}(e)}{2L^2 i} \sin^{-2} \left( \frac{z_S(e) - z_{v,j}}{L/\pi} \right) \right] \\
&+ \frac{\omega_0}{4\pi} \int_0^{L_S} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_S(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_S(e')}{L/(2\pi)} \right) \\
&\qquad \times \operatorname{Re} \left[ \partial_t z_{S,e}(e) \overline{z_{S,e}(e')} + z_{S,e}(e) \overline{\partial_t z_{S,e}(e')} \right] \mathrm{d}e' \\
&- \omega_0 \int_0^{L_S} \operatorname{Im} \left[ \frac{\partial_t z_S(e) - \partial_t z_S(e')}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] \operatorname{Re} \left[ z_{S,e}(e) \overline{z_{S,e}(e')} \right] \mathrm{d}e' \\
&- \frac{\omega_0}{4\pi} \int_0^{L_B} \ln \left( \cosh \operatorname{Im} \frac{z_S(e) - z_B(e')}{L/(2\pi)} - \cos \operatorname{Re} \frac{z_S(e) - z_B(e')}{L/(2\pi)} \right) \operatorname{Re} \left[ \partial_t z_{S,e}(e) \overline{z_{B,e}(e')} \right] \mathrm{d}e' \\
&+ \omega_0 \int_0^{L_B} \operatorname{Im} \left[ \frac{\partial_t z_S(e)}{2Li} \cot \left( \frac{z_S(e) - z_B(e')}{L/\pi} \right) \right] \operatorname{Re} \left[ z_{S,e}(e) \overline{z_{B,e}(e')} \right] \mathrm{d}e'
\end{align*}
$$

where we have used the relation on the cotangent (2.15).

Every integrals are classically defined, in particular the functions can be extended by continuity for $e = e'$ in the third integral by

$$
\operatorname{Re} \left[ \frac{\gamma_S(e) z_{S,ee}(e) - \gamma_{S,e}(e) z_{S,e}(e)}{2\pi i z_{S,e}^2(e)} \partial_t z_{S,e}(e) \right]
$$

and in the fourth one by

$$
\operatorname{Re} \left[ \frac{\gamma_S(e) \partial_t z_{S,ee}(e) - \gamma_{S,e}(e) \partial_t z_{S,e}(e)}{2\pi i z_{S,e}(e)} \right],
$$

which can be simplified by a part in the extension of the third integral. We can replace the
term $z_{S,ee}$ in the extension of the second integral by replacing $\frac{L_S}{N_S}$ Re $\left[ \frac{\gamma_S(i)z_{S,ee}(i)}{2\pi i z_{S,i}^2(e)} \partial_t z_{S,e}(i) \right]$
by

$$ \frac{L_S}{N_S} \sum_{j \neq i} \mathrm{Re} \left[ \frac{\gamma_S(i)}{z_{S,i}^2} \partial_t z_{S,e}(i) \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) z_{S,e}(j) \right]. $$

Unfortunately, it is more complicated to replace $\partial_t z_{S,e}$ and $\partial_t z_{S,ee}$ because differentiating the previous relation would introduce an additional $\sin^{-2}$ term. The first integral has to be replaced by

$$ \int_0^{L_S} \mathrm{Re} \left[ \frac{\partial_t \gamma_S(e') z_{S,e}(e)}{2Li} \cot \left( \frac{z_S(e) - z_S(e')}{L/\pi} \right) \right] de', $$

where the continuous function is extended for $e = e'$ by zero.

These computations provide the explicit expression for $G_{V,1}$. We can note that the many terms disappear when considering the single-fluids formulation, i.e. the $\alpha = 1$ and $A_{tw} = 1$ case.

For the expression of $G_{V,2}$, we get

$$ \begin{aligned} G_{V,2}(e) = & - \int_0^{L_S} \gamma_S(e') \mathrm{Im} \left[ \frac{\pi z_{B,e}(e) \partial_t z_S(e')}{2L^2 i} \sin^{-2} \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] de' \\ & - \sum_{j=1}^{N_v} \gamma_{v,j} \mathrm{Im} \left[ \pi \frac{z_{B,e}(e) \partial_t z_{v,j}}{2L^2 i} \sin^{-2} \left( \frac{z_B(e) - z_{v,j}}{L/\pi} \right) \right] \\ & - \frac{\omega_0}{4\pi} \int_0^{L_S} \ln \left( \cosh \mathrm{Im} \frac{z_B(e) - z_S(e')}{L/(2\pi)} - \cos \mathrm{Re} \frac{z_B(e) - z_S(e')}{L/(2\pi)} \right) \mathrm{Im} \left[ z_{B,e}(e) \overline{\partial_t z_{S,e}(e')} \right] de' \\ & - \omega_0 \int_0^{L_S} \mathrm{Im} \left[ \frac{\partial_t z_S(e')}{2Li} \cot \left( \frac{z_B(e) - z_S(e')}{L/\pi} \right) \right] \mathrm{Im} \left[ z_{B,e}(e) \overline{\partial_t z_{S,e}(e')} \right] de'. \end{aligned} $$

Concerning the numerical approximation, for $z_B(i) := z_B(e_{B,i})$, $z_S(i) := z_S(e_{S,i})$ and $\gamma_S(i) := \gamma_S(e_{S,i})$ given, we set $\tilde{z}_B(i) := (z_B(e_{B,i}) + z_B(e_{B,i+1}))/2$ for $i = 1, \dots, N_B - 1$ to construct the matrix $B_{B,N}$

$$ \begin{aligned} B_{B,N}(i,j) &= \frac{L_B}{N_B} \mathrm{Im} \left[ \frac{\tilde{z}_{B,e}(i)}{2Li} \cot \left( \frac{\tilde{z}_B(i) - z_B(j)}{L/\pi} \right) \right] \quad \forall (i,j) \in [1, N_B - 1] \times [1, N_B], \\ B_{B,N}(N_B, j) &= \frac{L_B}{N_B} \quad \forall j \in [1, N_B], \end{aligned} $$

The discretization of RHS$_{V0,B,N}$ and RHS$_{VB,N}$ are clear, replacing every $z_B(e)$ by $\tilde{z}_B(i)$, and where the last component is $-\gamma$. We deduce $\gamma_B = B_{B,N}^{-1}$ RHS$_{VB,N}$.

Next, we compute $\partial_t z_S$, $\hat{u}_F(z_S(e))$ and $\hat{u}_A(z_S(e))$ where we could extend in the integral on $\Gamma_S$ for $e' = e$ by $\frac{\gamma_S(e)z_{S,ee}(e) - \gamma_{S,e}(e)z_{S,e}(e)}{2\pi i z_{S,e}^2(e)}$, and we compute the derivative with respect to $e$.
We compute $A_{S,N}$ as

$$A_{S,N}(i,j) = A_{tw} \frac{L_S}{N_S} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_S(j)}{L/\pi} \right) z_{S,e}(i) \right] \quad \forall i \neq j \in [1, N_S] \times [1, N_S],$$

$$A_{S,N}(i,i) = \frac{1}{2} \quad \forall i \in [1, N_S],$$

next

$$C_{V,N}(i,j) = A_{tw} \frac{L_B}{N_B} \mathrm{Re} \left[ \frac{1}{2Li} \cot \left( \frac{z_S(i) - z_B(j)}{L/\pi} \right) z_{S,e}(i) \right] \quad \forall (i,j) \in [1, N_S] \times [1, N_B].$$

and finally

$$D_{V,N}(i,j) = \frac{L_S}{N_S} \mathrm{Im} \left[ \frac{1}{2Li} \cot \left( \frac{\tilde{z}_B(i) - z_S(j)}{L/\pi} \right) \tilde{z}_{B,e}(i) \right] \quad \forall (i,j) \in [1, N_B - 1] \times [1, N_S],$$

$$D_{V,N}(N_B, j) = 0 \quad \forall j \in [1, N_S].$$

## APPENDIX D. NOTATIONS

We summarize here the notations used in this article.

* For a given vector $\mathbf{a}_R = (a_{R,1}, a_{R,2})$, $a_{R,1}$ is the first component, whereas $a_{R,2}$ is the second component. We also introduce the complex notation: $a_R = a_{R,1} + ia_{R,2}$.

* For a function $e \mapsto f_R(e)$, we denote the derivative with respect to $e$ as $f_{R,e}$.

* We define three domains: $\mathcal{D} = \mathbb{T}_L \times \mathbb{R}$, $\mathcal{D}_F$ the fluid domain, $\mathcal{D}_A$ the air domain, $\mathcal{D}_B$ the domain below the bottom. $\Gamma_S$ is then the water-air free surface, and $\Gamma_B$ the bottom (see Fig. 1 in Section 1).

* Parametrization of the boundaries: the free surface $\Gamma_S$ (initially parameterised by arclength from left to right) $e \in [0, L_S] \mapsto z_S(t, e) = z_{S,1}(e, e) + iz_{S,2}(t, e)$; the bottom $\Gamma_B$ (by arclength from left to right) $e \in [0, L_B] \mapsto z_B(e) = z_{B,1}(e) + iz_{B,2}(e)$.

* Tangent vectors $\tau_S = \tau_{S,1} + i\tau_{S,2} = |z_{S,e}(e)|^{-1}z_{S,e}$ and $\tau_B = \tau_{B,1} + i\tau_{B,2} = |z_{B,e}(e)|^{-1}z_{B,e}$ are pointing to the right. The normal vector $n_S = n_{S,1} + in_{S,2} = -\tau_{S,2} + i\tau_{S,1} = i\tau_S$ is pointing out of the fluid domain, whereas $n_B = n_{B,1} + in_{B,2} = -\tau_{B,2} + i\tau_{B,1} = i\tau_B$ is pointing in the fluid domain.

* The jumps of the tangential components $(\gamma_S, \gamma_B)$ are defined in (2.12), whereas the jumps of the potentials $(\mu_S, \mu_B)$ are defined in (2.19).

* The mean current (or circulation) is $\gamma \in \mathbb{R}$ and the constant vorticity is $\omega_0 \in \mathbb{R}$. The non constant part of the vorticity inside the fluid is $\sum_{j=1}^{N_v} \gamma_{v,j} \delta_{z_{v,j}(t)}$.

* $\psi_F$ is the stream function associated to $\mathbf{u} = \nabla^\perp \psi_F$.

* $u_{\omega,\gamma}$ is defined as a vector satisfying (2.7) and (2.10), which allows to define $\mathbf{u}_R := \mathbf{u} - \mathbf{u}_{\omega,\gamma}$ as the gradient of a potential function (see (2.8)).

* For vector fields $\mathbf{u} = (u_1, u_2)$, we define $\hat{\mathbf{u}} = u_1 - iu_2$, $(u_1, u_2)^\perp = (-u_2, u_1)$, as well as the curl operator $\mathrm{curl} \, \mathbf{u} = \partial_1 u_2 - \partial_2 u_1$.
REFERENCES

[1] T. Alazard. Boundary observability of gravity water waves. *Ann. Inst. Henri Poincaré, Anal. Non Linéaire*, 35(3):751-779, 2018.

[2] T. Alazard, N. Burq, and C. Zuily. On the Cauchy problem for gravity water waves. *Invent. Math.*, 198(1):71-163, 2014.

[3] D. M. Ambrose, R. Camassa, J. L. Marzuola, R. M. McLaughlin, Q. Robinson, and J. Wilkening. Numerical algorithms for water waves with background flow over obstacles and topography. *Adv Comput Math*, 48, 2022.

[4] D. Arsenio, E. Dormy, and C. Lacave. The Vortex Method for Two-Dimensional Ideal Flows in Exterior Domains. *SIAM J. Math. Anal.*, 52(4):3881-3961, 2020.

[5] G. Baker. Generalized vortex methods for free-surface flows. In R. E. MEYER, editor, *Waves on Fluid Interfaces*, pages 53-81. Academic Press, 1983.

[6] G. R. Baker, D. I. Meiron, and S. A. Orszag. Generalized vortex methods for free-surface flow problems. *J. Fluid Mech.*, 123:477-501, 1982.

[7] G. R. Baker and C. Xie. Singularities in the complex physical plane for deep water waves. *J. Fluid Mech.*, 685:83-116, 2011.

[8] J. T. Beale, T. Y. Hou, and J. Lowengrub. Convergence of a boundary integral method for water waves. *SIAM J. Numer. Anal.*, 33(5):1797-1843, 1996.

[9] G. Beck and D. Lannes. Freely floating objects on a fluid governed by the Boussinesq equations. *Ann. Inst. Henri Poincaré, Anal. Non Linéaire*, 39(3):575-646, 2022.

[10] J. Beichman and S. Denisov. 2D Euler equation on the strip: stability of a rectangular patch. *Commun. Partial Differ. Equations*, 42(1):100-120, 2017.

[11] D. Bresch, D. Lannes, and G. Metivier. Waves interacting with a partially immersed obstacle in the boussinesq regime. *Anal. PDE*, 14(4):1085-1124, 2021.

[12] P. F. Byrd and M. D. Friedman. *Handbook of elliptic integrals for engineers and physicists*, volume 67. Springer, 2013.

[13] Y. Cao, W. W. Schultz, and R. F. Beck. Three-dimensional desingularized boundary integral methods for potential problems. *Int. J. Numer. Methods Fluids*, 12(8):785-803, 1991.

[14] A. Castro, D. Córdoba, C. Fefferman, F. Gancedo, and J. Gómez-Serrano. Finite time singularities for the free boundary incompressible Euler equations. *Ann. Math. (2)*, 178(3):1061-1134, 2013.

[15] W. Craig and C. Sulem. Numerical simulation of gravity waves. *J. Comput. Phys.*, 108(1):73-83, 1993.

[16] R. Dautray and J.-L. Lions. *Mathematical analysis and numerical methods for science and technology. Volume 4: Integral equations and numerical methods*. Berlin: Springer, 2nd printing edition, 2000.

[17] E. Fabes, M. Jodeit, and N. Riviere. Potential techniques for boundary value problems on $C^1$-domains. *Acta Math.*, 141:165-186, 1978.

[18] P. Germain, N. Masmoudi, and J. Shatah. Global solutions for the gravity water waves equation in dimension 3. *Ann. Math. (2)*, 175(2):691-754, 2012.

[19] J. Goodman, T. Y. Hou, and J. Lowengrub. Convergence of the point vortex method for the 2-D Euler equations. *Comm. Pure Appl. Math.*, 43(3):415-430, 1990.

[20] S. T. Grilli, J. Skourup, and I. Svendsen. An efficient boundary element method for nonlinear water waves. *Eng. Anal. Bound. Elem.*, 6(2):97-107, 1989.

[21] P. Guyenne and S. T. Grilli. Numerical study of three-dimensional overturning waves in shallow water. *J. Fluid Mech.*, 547:361-388, 2006.

[22] T. Iguchi and D. Lannes. Hyperbolic free boundary problems and applications to wave-structure interactions. *Indiana Univ. Math. J.*, 70(1):353-464, 2021.

[23] B. Jiang and Q. Bi. Classification of traveling wave solutions to the green-naghdi model. *Wave Motion*, 73:45-56, 2017.

[24] R. S. Johnson. *A modern introduction to the mathematical theory of water waves*. Camb. Texts Appl. Math. Cambridge: Cambridge University Press, 1997.

[25] O. Kellogg. *Foundations of potential theory*. Berlin-Heidelberg-New York: Springer-Verlag, 1967.

[26] D. Lannes. Well-posedness of the water-waves equations. *J. Am. Math. Soc.*, 18(3):605-654, 2005.
[27] D. Lannes. *The water waves problem. Mathematical analysis and asymptotics*, volume 188. Providence, RI: American Mathematical Society (AMS), 2013.

[28] G. Moon. Local well-posedness of the gravity-capillary water waves system in the presence of geometry and damping. *arXiv:2201.04713*, 2022.

[29] N. I. Muskhelishvili. *Singular integral equations*. Wolters-Noordhoff Publishing, Groningen, 1972. Boundary problems of functions theory and their applications to mathematical physics, Revised translation from the Russian, edited by J. R. M. Radok, Reprinted.

[30] Y. Pomeau, M. L. Berre, P. Guyenne, and S. Grilli. Wave-breaking and generic singularities of non-linear hyperbolic equations. *Nonlinearity*, 21(5):T61, 2008.

[31] Y. Pomeau and M. Le Berre. Topics in the theory of wave-breaking. In *Singularities in mechanics. Formation, propagation and microscopic description*, pages 125–162. Paris: Société Mathématique de France (SMF), 2012.

[32] A. Riquier and E. Dormy. A numerical study of the viscous breaking water waves problem and the limit of vanishing viscosity. In prep., 2023.

[33] Y.-M. Scolan. Some aspects of the flip-through phenomenon: A numerical study based on the desingularized technique. *J. Fluids Struct.*, 26(6):918–953, 2010.

[34] D. Scullen and E. Tuck. Nonlinear free-surface flow computations for submerged cylinders. *J. Ship Res.*, 39(3):185–193, 1995.

[35] E. Tuck. *Solution of Free-Surface Problems by Boundary and Desingularised Integral Equation Techniques*. World Scientific Publishing, 1992. Computational Techniques and Applications: CTAC97 B.J. Noye, M.D. Teubner and A.w. Gill (Editors).

[36] S. Wu. Well-posedness in Sobolev spaces of the full water wave problem in 2d. *Invent. Math.*, 130(1):39–72, 1997.

[37] V. E. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. *J. Applied Mech. Tech. Phys.*, 9(2):190–194, 1968.

(E. Dormy) DÉPARTEMENT DE MATHÉMATIQUES ET APPLICATIONS, UMR-8553, ÉCOLE NORMALE SUPÉRIEURE, CNRS, PSL UNIVERSITY, 75005 PARIS, FRANCE

Email address: Emmanuel.Dormy@ens.fr

(C. Lacave) UNIV. SAVOIE MONT BLANC, CNRS, LAMA, 73000 CHAMBÉRY, FRANCE.<fnref n="5" />

Email address: Christophe.Lacave@univ-smb.fr

<footnote n="5">⁵Previously at Univ. Grenoble Alpes, CNRS, IF, 38000 Grenoble, France</footnote>