# Air/water interfacial waves with a droplet at the tip of their crest

B C Nie, X Guan, J.-M Vanden-Broeck, F Dias

► To cite this version:

B C Nie, X Guan, J.-M Vanden-Broeck, F Dias. Air/water interfacial waves with a droplet at the tip of their crest. Physics of Fluids, 2023, 35 (1), pp.012101. 10.1063/5.0131944. hal-04497649

HAL Id: hal-04497649

https://hal.science/hal-04497649v1

Submitted on 10 Mar 2024

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 CC BY 4.0 - Attribution - International License
# Air/water interfacial waves with a droplet at the tip of their crest

Cite as: Phys. Fluids **35**, 012101 (2023); https://doi.org/10.1063/5.0131944

Submitted: 25 October 2022 • Accepted: 07 December 2022 • Accepted Manuscript Online: 09

December 2022 • Published Online: 03 January 2023

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

B. C. Nie (聂冰川),

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

X. Guan (关昕),

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

J.-M. Vanden-Broeck, et al.

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

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

View Online

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

Export Citation

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

CrossMark

## ARTICLES YOU MAY BE INTERESTED IN

*Some exact solutions that describe an ideal fluid flow with a free boundary*
*Physics of Fluids 35, 017104 (2023); https://doi.org/10.1063/5.0132100*

*Numerical simulation of secondary breakup of shear-thinning droplets*
*Physics of Fluids 35, 012103 (2023); https://doi.org/10.1063/5.0132939*

*Morphed inception of dynamic Leidenfrost regime in colloidal dispersion droplets*
*Physics of Fluids (2022); https://doi.org/10.1063/5.0131609*

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

## Physics of Fluids

### Special Topic: Paint and Coating Physics

**Submit Today!**

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

AIP
Publishing
# Air/water interfacial waves with a droplet at the tip of their crest

Cite as: Phys. Fluids **35**, 012101 (2023); doi: 10.1063/5.0131944

Submitted: 25 October 2022 · Accepted: 7 December 2022 ·

Published Online: 3 January 2023

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

View Online

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

Export Citation

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

CrossMark

B. C. Nie (聂冰川),<fnref n="1" />,<fnref n="2" /> ID X. Guan (关昕),<fnref n="3" /> ID J.-M. Vanden-Broeck,<fnref n="3" /> ID and F. Dias<fnref n="2" />,<fnref n="4" />,α) ID

## AFFILIATIONS

## ABSTRACT

In nature, it is common to observe water wave crests with a droplet at their tip. This fascinating configuration remains unexplained from the physical point of view. The present study explores such a unique local configuration numerically. Solitary waves that propagate at the interface between two layers of irrotational fluids are considered. Extending the work of Guan et al. ["A local model for the limiting configuration of interfacial solitary waves," J. Fluid Mech. 921, A9 (2021)], the density ratio has been decreased to a very small value equal to 0.001, which is close to the air/water density ratio at sea level (0.0013). A highly accurate solution for the limiting configuration of solitary waves is obtained by solving the irrotational Euler equations using the boundary integral method and Newton's iterations. It is confirmed that the limiting configuration consisting of a droplet sitting on a crest with a 120° angle exists for very small density ratios. This limiting configuration obviously does not exist for surface waves with a void on the top, thus stressing the crucial role played by the air. The droplet is stationary in a frame of reference moving with the wave and experiences intense shear at its tip. From the energy point of view, the formation of a crest with a droplet is accompanied by a remarkable drop of kinetic and potential energies of water in the vicinity of the crest. Furthermore, we present a simple set of scaling relations for the fall of the droplet.

© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0131944

## I. INTRODUCTION

Water waves with a droplet on their crest can be observed on the open ocean, on the seashore, and on lakes. Figure 1 shows some examples observed on the seashore in south Dublin, Ireland, and one may also refer to the field observation in Ref. 1 (Fig. 5). This common, yet unexplained, phenomenon may be related to wave breaking as these droplets quickly bend and splash on the sea surface, thus generating foam and spume. At the end of their experimental study, Villermaux and Pomeau<fnref n="1" /> commented that the formation of droplets on the crest may be relevant to the superfall of wave crests. Superfall means here that the fluid can experience accelerations that are larger than the acceleration due to gravity only. However, their experiments were carried out with a different model than that of water waves, i.e., free fall of water filling a vertical tube with expanding cross-section and open ends. They did observe a concentrated "nipple" on the top of the essentially flat base solution and attributed it to the instability of the interface due to a pressure gradient during the superfall.

In fact, the droplet-type crest could be one of the limiting configurations of water waves. The investigation of the limiting configuration of wave crests goes back to the conjecture made by Stokes one and a half centuries ago about the crest profile of periodic traveling waves when they reach their maximal amplitude. They are known as the Stokes highest waves. After many years of discussion from asymptotic, numerical, and analytical perspectives,<fnref n="2" />⁻<fnref n="4" /> it is now widely accepted that the limiting configuration of the Stokes highest wave is a stagnation point at the crest with an enclosed angle of 120°. At the same time, the limiting configuration of solitary waves has also been investigated. It has been shown that extreme solitary waves have the same limiting configuration as the periodic waves.⁵⁻⁸ Indeed, the limiting configuration at the very tip is localized and independent of wavelength and water depth.<fnref n="2" /> Comparing this limiting configuration obtained with a surface wave model to the droplet-type crest, it looks like the surface wave model can hardly reproduce and explain the formation of a droplet on the crest. One possible reason is that the

<footnote n="1">¹School of Physical Science and Engineering, Beijing Jiaotong University, Beijing 100044, China</footnote>

<footnote n="2">²School of Mathematics and Statistics, University College Dublin, Dublin, Ireland</footnote>

<footnote n="3">³Department of Mathematics, University College London, London WC1E 6BT, United Kingdom</footnote>

<footnote n="4">⁴Centre Borelli, ENS Paris-Saclay, Gif-sur-Yvette 91190, France</footnote>

<footnote>a) Author to whom correspondence should be addressed: frederic.dias@ucd.ie</footnote>
<figure><img src="image_14.png" /><figcaption>FIG. 1. Wave crests with a droplet at their tip observed on the seashore in south Dublin, Ireland: (a) an example of a dense distribution and (b)–(g) view of the typical local forms from different angles.</figcaption></figure>

influence of air which has been neglected so far turns out to be more
important than originally thought.

In this paper, we take air into account and investigate the limiting
configuration of interfacial waves between two layers of homogeneous
fluids. The aforementioned sharp crest of 120° is no longer possible.
Otherwise, the fluid velocity in the upper layer near the crest would
become infinite.⁹ A wave profile with a vertical tangent was initially
regarded as the configuration of the highest-possible interfacial
wave.¹⁰ Later, the existence of overhanging configurations was found
numerically.¹¹,¹² Meiron and Saffman⁷ asserted that these overhanging
solutions would continue to exist until the interface becomes self-
intersecting. Finally, a stagnant fluid “bubble” on the top of a 120°
angle was proposed by Grimshaw and Pullin¹³ as the limiting configu-
ration of interfacial waves. Recently, a highly accurate numerical inves-
tigation of the limiting configuration of interfacial waves was
conducted (see Ref. 14, Fig. 7), in which the smallest density ratio used
was 0.01. This configuration being localized, it is also valid for solitary
interfacial waves.¹⁵,¹⁶ In fact, the limiting configuration of a stagnant
fluid bubble on a wave crest is highly similar to the droplet-type crest.
This similarity led us to wonder about their connections.

The density ratio ρ could be the most crucial number for the
localized configuration of the crest under the precondition that the
depth allows the wave to develop freely. However, the aforementioned
investigations are restricted to moderate values of ρ. Although

interfacial waves with ρ = 10⁻⁴ were examined in Ref. 17, the over-
hanging profile with the smallest value of ρ is only given for ρ = 0.15.
Recently, Guan et al.¹⁸ revisited the limiting configuration of solitary
waves and decreased the value of ρ to 0.01. To overcome the rapid
shrinking of the local limiting configuration when decreasing ρ even
further, a simplified model was proposed. In that study, it was found
that the formation of the limiting configuration is always accompanied
by a decrease in the amplitude. It seems to be consistent with the con-
cept that nipples are caused by the superfall of the crest, see Ref. 1.
Considering that the actual density ratio of air to water is of the order
of O(10⁻³), that is at least one order of magnitude smaller than the
values used in the existing studies. For such a low density ratio, the
existence of a limiting configuration with a stagnant fluid bubble on
the crest is still an open question. If this limiting configuration does
exist, the actual solution has not been presented yet, and the possible
relationship with the droplet-type crests has not been discussed yet.

In the present study, solitary interfacial waves are investigated.
We follow Ref. 18 and decrease the value of ρ to the air–water density
ratio (\~10⁻³). It turns out that the limiting configuration with a stag-
nant droplet on the wave crest does exist for such a low density ratio.
The actual configuration, the velocity field and the shear along the
interface are presented. The variation in the kinetic and potential ener-
gies of the fluid in the crest region along the solution branch is dis-
cussed. It is shown that the decrease in the wave amplitude is not a
necessary condition for reaching the limiting configuration. Instead, the remarkable drop of mechanical energy inside the crest is. Furthermore, a simple set of scaling relations is introduced for estimating the splash of the main droplet. Our results may be helpful for understanding the common, yet unexplained, droplet-type crest phenomenon.

## II. FORMULATION OF THE PROBLEM

We consider a two-dimensional air layer of density $ρ_2$ and thickness $h_2$ over a water layer of density $ρ_1$ and thickness $h_1$, as illustrated in Fig. 2. The top and bottom boundary conditions are solid horizontal walls. The stratification is overall stable as $ρ_1 > ρ_2$. An interfacial solitary wave travels at speed c along the x-axis from the left ($x → -∞$) to the right ($x → ∞$), and it is symmetric with respect to the vertical y-axis. The system of coordinates moves with the wave, and the level $y=0$ is the undisturbed level of the interface. Gravity acts in the negative y direction. The quantities $h_1$, c, and $ρ_1$ are used as length, velocity, and density units, respectively, to nondimensionalize the system.

Under the assumption of irrotational flow, the motion of water and air can be described by the Laplace equation

$$ \frac{\partial^2 \phi_j}{\partial x^2} + \frac{\partial^2 \phi_j}{\partial y^2} = 0 \quad \text{in } \Omega_j, \quad j = 1, 2. \qquad (1) $$

$\phi_j$ is the nondimensional velocity potential, and $\Omega_j$ is the fluid domain shown in Fig. 2. Here and hereafter, the subscripts j=1 and 2 correspond to the water layer and the air layer, respectively.

The two kinematic and the dynamic boundary conditions on the interface read

$$ \frac{\partial \phi_j}{\partial y} - \frac{\partial \phi_j}{\partial x} \frac{\partial \eta}{\partial x} = 0, \quad j = 1, 2 \qquad (2) $$

and

$$ \rho |\nabla \phi_2|^2 - |\nabla \phi_1|^2 + 2(\rho - 1) \frac{\eta}{F^2} = \rho - 1. \qquad (3) $$

Here, $ρ = ρ_2/ρ_1$ represents the air/water density ratio, $η$ is the nondimensional elevation of the interface, $F = c/\sqrt{gh_1}$ is the Froude number, and $g$ is the acceleration due to gravity.

<figure><img src="image_15.png" /><figcaption>FIG. 2. Schematic of the crest of a solitary wave on the air-water interface. The wave $y = η(x)$ propagates at a constant velocity c. It lies between two layers of fluids with density $ρ_1$ and depth $h_1$, i=1, 2, with $ρ_1 > ρ_2$. The crest domain $Ω'$ is a square of length $l$. The wave shown is a numerical solution $y = η(x)$ with $ρ = 0.001, h = 80$.</figcaption></figure>

The two fluid layers are bounded by the horizontal upper and bottom walls. The corresponding boundary conditions can be expressed as

$$ \begin{aligned} \frac{\partial \phi_1}{\partial y} &= 0 \quad \text{at } y = -1, \\ \frac{\partial \phi_2}{\partial y} &= 0 \quad \text{at } y = h. \end{aligned} \qquad (4) $$

The depth ratio $h$ is defined as $h_2/h_1$.

In the frame of reference moving with the wave speed, the boundary conditions at infinity for the solitary wave are

$$ \frac{\partial \phi_j}{\partial x} \to -1, \quad \frac{\partial \phi_j}{\partial y} \to 0, \quad \eta \to 0 \quad \text{at } x \to \pm\infty, \quad j = 1, 2. \qquad (5) $$

The solitary waves are characterized by three dimensionless parameters: $h$, $ρ$, and $F$. The impact of $h$ is insignificant when it is large so that the wave can develop freely. According to Kataoka (Ref. 17, Table I), the threshold value is of the order of $h > 5$ for small $ρ$. In this study, $h$ is fixed at 80. This implies that when $ρ$ will be further fixed at $10^{-3}$, a value close to the air/water density ratio, there will be only one governing parameter: the Froude number $F$.

## III. NUMERICAL METHOD

The numerical method closely follows the numerical method used in Ref. 18. Some details are repeated for the sake of completeness, and the differences are emphasized. The interface is parameterized by the arc length $s$ as $x = X(s)$ and $y = η(s)$. It gives the arc length equation

$$ X'(s)^2 + \eta'(s)^2 = 1. \qquad (6) $$

Let $s=0$ at $x=0$. There are two conditions for $X$ and $η$

$$ X(0) = 0, \quad η(0) = A. \qquad (7) $$

Here, $A$ is the nondimensional wave amplitude.

Following Refs. 18 and 19, the Laplace equation (1) can be reformulated using the Cauchy integral formula. Let $u$ and $v$ denote the horizontal and vertical velocities, respectively. Introducing $Ψ_j(z) = u_j(z) - iv_j(z) + 1$, which is an analytic function of $z = x + iy$, the Cauchy integral formula gives

$$ \Psi_j(z_0) = \frac{1}{\pi i} \oint_{\partial \Omega_j} \frac{\Psi_j(z)}{z - z_0} dz, \quad j = 1, 2, \qquad (8) $$

for $z_0$ on the boundary of the integration domain $Ω_j$.

We use the Schwarz reflection principle to satisfy the boundary conditions on the upper and bottom walls (4). The integration domain in (8) for water can be $Ω_1$ plus its reflection about the bottom wall, which is bounded by $η$ and $η_{R,1}$. Similarly, for the air, its integration domain $Ω_2$ is bounded by $η$ and $η_{R,2}$. $η_{R,j}$ are expressed as

$$ \eta_{R,j} = \begin{cases} -2 - \eta & \text{when } j = 1, \\ 2h - \eta & \text{when } j = 2. \end{cases} \qquad (9) $$

Based on the arc length parameterization, the integration in (8) carried out in the original plus reflected domains can be rewritten as
$$ \Psi_j(s_0) = \frac{1}{\pi i} \oint_{\eta+\eta_{R,j}} \frac{\Psi_j(s)[X' + iY']}{[X(s) - X(s_0)] + i[Y(s) - Y(s_0)]} ds, \quad j = 1, 2. \quad (10) $$

According to (5), $\Psi_j \to 0$ when $x \to \infty$. Thus, the integral above can be calculated in $s \in (-L, L)$ with large $L$. Using the symmetry of the interface with respect to the y-axis and taking the real part, (10) can be finally simplified as

$$ u_j(s_0) + 1 = + \int_0^L \frac{[(u_1(s) + 1)X'(s) + v_1(s)\eta'(s)][\eta(s) - \eta(s_0)] - \eta'(s)X_-}{X_-^2 + [\eta(s) - \eta(s_0)]^2} ds + \int_0^L \frac{[(u_1(s) + 1)X'(s) + v_1(s)\eta'(s)][\eta(s) - \eta(s_0)] - \eta'(s)X_+}{X_+^2 + [\eta(s) - \eta(s_0)]^2} ds \\ - \int_0^L \frac{[(u_1(s) + 1)X'(s) + v_1(s)\eta'(s)][\eta_{R,j}(s) - \eta(s_0)] + \eta'(s)X_-}{X_-^2 + [\eta_{R,j}(s) - \eta(s_0)]^2} ds - \int_0^L \frac{[(u_1(s) + 1)X'(s) + v_1(s)\eta'(s)][\eta_{R,j}(s) - \eta(s_0)] + \eta'(s)X_+}{X_+^2 + [\eta_{R,j}(s) - \eta(s_0)]^2} ds. \quad (11) $$

Here, $X_\pm = X(s) \pm X(s_0)$ and $j = 1, 2$.

Meanwhile, the kinematic boundary condition (2) can also be written in the parameterized form

$$ v_j(s)X'(s) - u_j(s)\eta'(s) = 0, \quad j = 1, 2. \quad (12) $$

The integration domain $s \in (0, L)$ is divided into $N-1$ elements by seeds $s_k$, where $k = 1, 2, 3, \dots, N$. The length of the kth element $ds_k$ is equal to $s_{k+1} - s_k$, and the middle point of each element $\sigma_k$ is equal to $(s_{k+1} + s_k)/2$, where $k = 1, 2, 3, \dots, N-1$. In Ref. 18, uniform meshes were used. However, as the size of the droplet-type crest decreases significantly with the decrease in $\rho$, nonuniform meshes are adopted in the integral domain in the present paper. More specifically, uniform finer meshes for the crest region, uniform larger meshes in the tail, and meshes with fixed stretching ratio $ds_k/ds_{k-1} = 1.004$ in between the crest and tail are used.

The midpoint rule is applied for the integral equations (11) to avoid the singularities. It gives $2N-2$ algebraic equations. The boundary conditions of the interfaces (3), (6), and (12) are discretized on the seeds, which gives $2N$, $N$, and $N$ algebraic equations, respectively. Using the boundary condition (7), the number of unknown variables can be reduced to $6N$, which are $x'(s_k)$, $\eta'(s_k)$, $u_j(s_k)$, $v_j(s_k)$, $k = 1, 2, 3, \dots, N$. Together with the conditions $\eta'(0) = 0$ and $u_1(L) = -1$, we are ready to search the branch of steady solutions in the $F-A$ plane.

Lagrange second-order interpolation is adopted for the derivatives. Iterations are carried out by Newton's method with a maximum residual error less than $10^{-9}$. When the first turning point on the

branches appears, Lagrange third-order interpolation and a maximum residual error threshold of $10^{-11}$ are adopted.

The solution branches are explored by the following procedures: (i) a Gaussian profile with infinitesimal amplitude is produced as an initial guess of the solitary wave, and then, the first solution, a point on the branch, can be obtained by iteration; (ii) the former solution is used as the initial guess for the next step with a slight change of $A$ (or $F$); thus, the next new solution can be obtained by iteration; (iii) repeat step (ii) until the droplet-type configuration develops; and (iv) the solution branch in the $F-A$ plane is then built by collecting all the solutions.

Figure 3 shows that the mesh size at the crest $ds_1$ has a significant impact on the crest profile, whereas the domain size $L$ mainly affects the profile of the tail of the wave and the wave amplitude. Using 1900 elements with $ds_1 = 1.0 \times 10^{-3}$ and $L = 9.2$ was found to be a good balance between computational cost and precision.

# IV. RESULTS AND DISCUSSION

## A. Limiting configuration

Figure 4 presents the branch of solitary wave solutions in the $F-A$ plane. Eight typical solutions marked by the red circles along the branch are selected. They are S1 (1.0957, 0.2042), S2 (1.1830, 0.4125), S3 (1.2558, 0.6166), S4 (1.2974, 0.8095), S5 (1.2949, 0.8320), S6 (1.2916, 0.8530), S7 (1.2922, 0.8590), and S8 (1.2929, 0.8608). Halves of the symmetric crests are plotted in Fig. 5. It is shown that the angle of the crest decreases and converges to 120° along the branch. Figure 5(b)

<figure><img src="image_16.png" /><figcaption>FIG. 3. Sensitivity test of the mesh size $ds$ and the integral domain $L$ on the profile of the interface at the crest and in the tail: (a) wave crest and (b) wave tail.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 4. Solution branch in the F-A plane for $\rho = 0.001$ and $h = 80$. The points labeled S1-S8 correspond to typical solutions selected to observe the wave profile and the shear at the interface.</figcaption></figure>

confirms that the droplet-type crest does exist for the small density ratio $\rho = 0.001$ and that it develops after the second turning point S6.

Figure 4 also shows that the wave amplitude increases monotonically along the solution branch. However, for the solution branches with larger values of the density ratio $\rho$ (see Fig. 14 for $\rho = 0.1$ and illustrations for other density ratios in Ref. 18), the amplitudes always steepen along the solution branch to reach the highest wave and then decrease to form the droplet-type crest. Moreover, Villermaux and Pomeau¹ conjectured that the nipples on the wave crest may be caused by its superfall. It may give an illusion that the formation of the droplet-type crest is always accompanied by the decrease in A. However, the case $\rho = 0.001$ considered in the present study indicates that a decrease in A is not a necessary condition for the formation of the limiting configuration.

To extend the solution branch past S8, much finer meshes are needed as the length of the neck region of the droplet crest will decrease further. In Ref. 18, a simplified model was introduced: the crest was replaced by two straight solid walls intersecting at the origin with an angle of 120° and a closed fluid bubble with a flat bottom on the top of the angle. It was found that the profiles for different values of the density ratio $\rho$ are geometrically similar. The limiting

configuration obtained with the local model is also presented for $\rho = 0.001$ in Fig. 5(b). Comparing the S8 profile with the local model, it is seen that the neck region may shrink further, leading to self-intersection after S8 and the local model may still be a good representation of the true solution.

Next, we explore the velocity fields. The velocity components in the moving frame at any location $z_i$ in $\Omega_{j, j=1, 2}$, are calculated again using the Cauchy integral formula, yielding

$$ 
\begin{aligned}
u_j(z_i) &= \operatorname{Re} \left[ \frac{1}{2\pi i} \oint_{\eta+\eta_{R,j}} \frac{\Psi_j(z)}{z-z_i} dz \right] - 1, \\
v_j(z_i) &= -\operatorname{Im} \left[ \frac{1}{2\pi i} \oint_{\eta+\eta_{R,j}} \frac{\Psi_j(z)}{z-z_i} dz \right].
\end{aligned} 
\quad (13) $$

Let us take S4 and S8 as solutions that are representative of the classical and droplet-type crests, respectively. Their velocity fields with the wave and interfacial shear velocity are plotted in Fig. 6. The components of the interfacial shear velocity ($\Delta u, \Delta v$) are ($u_1 - u_2, v_1 - v_2$). Figures 6(a) and 6(c) show that the most intense shear occurs at the tips because of the significant slowdown of the water (in the moving frame) and speedup of air around the crest, and the shear velocity can be more than twice the wave speed. For the droplet-type situation, Fig. 6(b) indicates that the droplet at the wave tip is almost stationary. The most severe shear still occurs at the crest. However, it has been further strengthened to be more than five times larger than the wave speed, see Fig. 6(d). It is interesting that there exists a segment on the interface without shear, denoted by $S_b$, which corresponds to the segment labeled by the blue line at the bottom of the droplet, see Fig. 6(b). In fact, the velocity components at segment $S_b$ vanish as well.

The ratio between the speed of the water particle at the tip and the wave speed, $u_{1f} + 1$, along the solution branch is plotted in Fig. 7(a). It can be seen that the speed of water at the tip approaches the wave speed with an obvious speedup after point S4 along the solution branch. According to the kinematic wave breaking criterion (e.g, Ref. 20), it suggests that breaking won't occur for the droplet-type crest if 1.0 is regarded as the threshold ratio between the particle speed at the tip and the wave speed. Note that the value of 1.0 was revisited by Barthelemy et al.²¹ and brought down to 0.85, in which case breaking would occur. Figure 7(b) shows the square of the shear velocity at the tip, along the solution branch. It seems that there is a simple exponential relation

<figure><img src="image_18.png" /><figcaption>FIG. 5. Shape of the half symmetric wave crest of eight selected typical solutions S1-S8 (see Fig. 4). A is the wave amplitude, $\rho = 0.001$, and $h = 80$: (a) waves S1, S2, S3, and S4; (b) waves S5, S6, S7, S8, and solution obtained with the local model in Ref. 18 with $\rho = 0.001$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 6. Velocity fields in the frame moving with the wave and shear velocity at the interface, for $\rho = 0.001$ and $h=80$: (a) and (c) correspond to solution S4 with a normal crest, and (b) and (d) correspond to solution S8 with a droplet-type crest. $S_b$ denotes the segment on the interface without shear.</figcaption></figure>

$$ (u_{2t} - u_{1t})^2 \sim A^{\alpha}, \qquad (14) $$

with $\alpha = 2.1$ for small and moderate amplitude and $\alpha = 42.5$ for extreme high amplitude.

## B. Crest mass and energies

The kinetic and potential energies in a domain $\Omega'$ around the crest, called the “crest domain” hereafter, are of interest as the limiting configuration is localized. The kinetic energy $E_{T,j}$, the potential energy $E_{P,j}$, and the mass $M_j$ of water and air in that domain are calculated as

$$ E_{T,j} = \frac{1}{2} \int \int_{\Omega'_j} \rho_j [(u_j + 1)^2 + v_j^2] F^2 dxdy, \qquad (15) $$

$$ E_{P,j} = \iint_{\Omega'_j} \rho_j y dxdy, \qquad (16) $$

$$ M_j = \iint_{\Omega'_j} \rho_j dxdy, \qquad (17) $$

where $j = 1, 2$. It is noteworthy that $E_{T,j}$ is calculated in the still frame of reference, and all the energies are nondimensionalized by $\rho_1 g h_1^3$.

The mass and energy can be influenced by the shape and size of the crest domain. In practice, the crest domain is a square domain of length $l$, and the wave tip is located at the height $5/6l$ and the horizontal center as marked by the domain $\Omega'$ in Fig. 2. Five different values of $l$, namely, 0.06 [as shown by the domain in Fig. 6(b)], 0.12, 0.24, 0.48, and 0.96, are examined. The results indicate that the variation

<figure><img src="image_20.png" /><figcaption>FIG. 7. Ratio between the speed of the water particle at the tip of the crest and wave speed (a) and shear velocity at the very tip of the wave (b). $u_t = u_j$ ($s=0$), $j=1, 2$, $\rho = 0.001$ and $h=80$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 8. Mass of water in the crest domain along the solution branch. $\rho = 0.001$ and $h = 80$. The length $l$ of the crest domain is set equal to 0.06. The crest corresponding to solution S8 is shown in Fig. 5(b).</figcaption></figure>

trends of $M_j$, $E_{T,j}$, and $E_{P,j}$ along the solution branch are insensitive to $l$, but not the locations of their extrema. The results corresponding to the crest domain with $l = 0.06$ are presented first.

With increasing wave amplitude along the solution branch shown in Fig. 4, the wave crest steepens from an initially flat interface and finally develops into a droplet-type crest. This means that part of the water in the crest domain is replaced by air, resulting in the decrease in $M_1$. Figure 8 shows that up to 38% of water has been lost along the solution branch. The sharp decrease starts at S4, the first turning point in the solution branch.

In general, the kinetic and potential energies of water and air are divided into two stages along the solution branch, i.e., a gentle variation stage followed by a sharp rise/drop stage, see Figs. 9 and 10. However, the turning points are not synchronized exactly. In fact, the size of the crest domain $l$ can also influence the location of the turning point slightly. Figure 9 indicates that $E_{T,1}$ and $E_{P,1}$ drop significantly during the formation of the droplet-type crest, although the wave amplitude does not decrease.

With the increase in the velocity of water at the tip, the kinetic energy of water in the crest domain first increases. The sharp drop of $E_{T,1}$ in Fig. 9(a) between S6 and S8 is mainly due to the remarkable decrease in the water mass in the crest domain, see Fig. 8. This implies that the variation in the kinetic energy of water is a result of the competition between an increase in the velocity and a decrease in the mass

in the crest domain. The sharp drop of $E_{P,1}$ in Fig. 9(b) is also due to the decrease in $M_1$.

The monotonic increase in $E_{T,2}$ and $E_{P,2}$ is straightforward, because the velocity of air at the tip keeps increasing up to more than five times larger than the wave speed [see Fig. 6(d)], and up to 38% of water is monotonically replaced by air (see Fig. 8) along the solution branch. The insignificant difference of $E_{T,2}$ between S7 and S8 is due to the flow of air by the side of the neck, i.e., the region between the profiles of S7 and S8 [see Fig. 5(b)] is characterized by small velocities.

The total energy of water and air, $E = \Sigma E_{T,j} + \Sigma E_{P,j}$, in the crest domain (“crest energy” hereafter) is plotted vs $F$ in Fig. 11. The maximum value of the crest energy as a function of $F$ is marked with a red cross. Figures 11(a)–11(d) show that the location of the maximum of the crest energy is sensitive to the size of the crest domain and shifts from S4 to S6 when $l$ decreases. This is because the energy of the smaller crest domain tends to reflect more the energy variation of the very tip. The droplet-typed crest formed at the very tip after S6 is accompanied by a significant energy drop. Overall, the formation of a droplet-type crest won’t occur before that maximum.

Solitary waves can be regarded as waves with infinite wavelength. Therefore, the wavelength of any disturbance will be shorter than the solitary wave. The stability of the basic wave to disturbances with a smaller wavelength is called superharmonic instability.²² The exchange of superharmonic instability occurs at stationary values of the total wave energy relative to the wave speed, as first found numerically by Tanaka²³ for periodic gravity waves in deep water. It was then analytically proved²⁴ and extended to surface solitary waves.²⁵ However, the situation is subtle for waves in water of finite depth. Indeed, Kataoka²⁶ showed that the exchange of stability for surface gravity waves occurs when the total wave energy becomes stationary as a function of the wave speed for fixed mean surface height and not for fixed Bernoulli’s constant. Kataoka¹⁷ showed that the energy criterion is also valid for interfacial waves. The application of the criterion for exchange of superharmonic instability to the branch of solutions computed in the present work is left for future work that will be fully dedicated to the stability of the present solitary waves.

We now present the results on the total energy in the whole computation domain as they are interesting in themselves, even without mentioning stability. This energy, which we call the “global energy” hereafter, is calculated as

$$E_G = \Sigma E_{GT,j} + \Sigma E_{GP,j} \quad \text{in } \Omega_j, j = 1, 2. \quad (18)$$

$E_{GT,j}$ and $E_{GP,j}$ are the global kinetic energy in the still frame of reference and global potential energy, respectively. Using the kinematic

<figure><img src="image_22.png" /><figcaption>FIG. 9. Energies of water in the crest domain: (a) kinetic energy measured in the still frame and (b) potential energy, the reference being the undisturbed free surface. The parameters are the same as in Fig. 7.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. 10. Energies of air in the crest domain: (a) kinetic energy measured in the still frame and (b) potential energy, the reference being the undisturbed free surface. The parameters are the same as in Fig. 7.</figcaption></figure>

boundary condition and (6), and nondimensionalizing by $\rho_1 gh_1^3$, the energies are expressed as

$$ 
\begin{aligned}
E_{GT,j} &= (-1)^{j-1} \frac{\rho_j}{2} F^2 \int_{-L}^{L} (\phi_j + x) \left. \frac{\partial(\phi_j + x)}{\partial n} \right|_{y=\eta} ds \\
&= (-1)^j \rho_j F^2 \int_{0}^{L} \left( \int_{0}^{s} (u_j(s)X'(s) + V_j(s)\eta'(s))ds + X(s) \right) \eta' ds
\end{aligned}
\quad (19) $$

and

$$ E_{GP,j} = (-1)^{j-1} \rho_j \int_{0}^{L} \eta(s)^2 X'(s) ds. \quad (20) $$

Here, $n$ is the normal to the interface.

The global energy along the branch with fixed Bernoulli's constant is presented in Fig. 12. The global energy is insensitive to $L$ as

long as $L$ is large enough. Indeed, $u|_{x=L} = v|_{x=L} \sim 0$, so there is no contribution to the kinetic energy. As for the potential energy, it is on the undisturbed free surface, and the constant term has been dropped during integration. The first stationary point ($dE_G/dF = 0$) occurs before S4, and the second stationary point occurs around S6.

Comparing Fig. 12 with Fig. 11, we conclude that the energy will first reach its maximum globally and then locally at the tip along the solution branch.

### C. Fall of the droplet

So far, we considered only stationary solutions in the frame of reference moving with the wave. In this section, we imagine that we consider time-dependent solutions. As time increases, the droplet will fall on the front side of the wave under gravity, considering the weak support from the neck, and finally, it will splash on the interface forming foam resulting in intense energy transfer and

<figure><img src="image_24.png" /><figcaption>FIG. 11. Maximum values (red crosses) of the total wave energy in the crest domain with different sizes of the crest domain. $l = (a) 0.06$, $(b) 0.12$, and $(c) 0.24$. $(d)$ is the location of the maximum value of $E$ on the branch when varying $l$. $\rho = 0.001$ and $h = 80$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIG. 12. Global wave energy $E_G$ as a function of the Froude number $F$ in the whole computation domain. $\rho = 0.001$ and $h = 80$.</figcaption></figure>

dissipation. At the same time, the wave will experience stretching as shown by the red profile in the sketch of Fig. 13. The stretched water body may turn into a thin film. This thin film can encounter significant aerodynamic effects and air entrainment and may further break into some micro-droplets released to the air.

There exist two different scenarios: one is that the droplet can fall a distance approximately equal to $A$, designated as *uninterrupted fall*, and the other one is that the droplet splashes on the slope profile of the wave itself with a fall that is much smaller than $A$, designated as *interrupted fall*. During the short period of the fall of the droplet, we introduce the following assumptions: (i) the limiting configuration of the wave consists of a droplet on the top, a main water body of the wave which is confined by two flat interfaces intersecting with an angle of 120°, and a neck region in between connecting them, as illustrated by the black dash line in Fig. 13; (ii) the droplet falls freely under gravity, and the aerodynamic force is much smaller than the gravity force; and (iii) the main water body of the wave propagates at a speed much slower than the wave speed with insignificant deformation. In fact, the last two assumptions imply that intense stretching would occur at the

<figure><img src="image_26.png" /><figcaption>FIG. 13. Sketch of the time-dependent evolution of the droplet-type crest.</figcaption></figure>

neck region. Based on the assumption of free fall, the timescale for the main droplet to fall a distance $s_v$ is $(2s_vh_1/g)^{1/2}$, where $s_v$ is normalized by $h_1$ and $0 < s_v < A$. If $c(2s_vh_1/g)^{1/2}$ is always larger than $s_vh_1 \tan(\pi/3)$, it goes to the uninterrupted fall. Otherwise, it goes to the interrupted fall. That gives a criterion identifying these two scenarios, which reads

$$A \lesssim \frac{2}{3} F^2. \qquad (21)$$

The criterion (21) has been plotted in the $F-A$ plane with the solution branches of $\rho = 0.001$ and $0.1$, as shown in Fig. 14. It suggests that the main droplet will experience an uninterrupted fall for the case of $\rho = 0.001$, whereas an interrupted fall of the main droplet can occur for a larger density ratio as indicated by the solution branch of $\rho = 0.1$.

Using the fact that the droplet falls a distance $A$ without interruption, the length scale of the droplet and its impacting velocity and angle can be further estimated. As discussed in Fig. 6(d), we have $u_{1t} = v_{1t} = v_{2t} = 0$ at the top of the droplet ($S_a$) and $u_1 = v_1 = u_2 = v_2 = 0$ at the bottom ($S_b$, see Fig. 13). Accordingly, the dynamic boundary condition (3) can be simplified, which gives the vertical coordinates of $S_a$ and $S_b$ as $1/2F^2 + 1/2\rho u_{2t}^2 F^2/(1 - \rho)$ and $1/2F^2$, respectively. Their difference gives the following dimensional length scale $l_s$ for the droplet:

$$l_s \sim \frac{\rho}{2(1 - \rho)} u_{2t}^2 F^2 h_1. \qquad (22)$$

This relationship for the length scale is valid as long as the droplet-type solution exists, whatever the value of the density ratio is. However, measuring the speed $u_{2t}$ is usually not easy. We may estimate $u_{2t}$ using the exponential relation in (14), $u_{2t} \sim A^{\alpha}$, with $\alpha = 42.5$, according to Fig. 7(b).

Based on the free fall timescale $(2Ah_1/g)^{1/2}$ and on the wave speed, we can also obtain the scales of the dimensional impact velocity $V_s$ and of the angle $\theta_s$. They read

<figure><img src="image_27.png" /><figcaption>FIG. 14. Regions for the fall of the main droplet. Uninterrupted fall: $A \lesssim 2/3F^2$; interrupted fall: $A \gtrsim 2/3F^2$.</figcaption></figure>
$$V_s \sim \sqrt{2gAh_1 + c^2}, \\ \theta_s \sim \tan^{-1}\left(\sqrt{gAh_1/c}\right). \qquad (23)$$

The combination of (22) and (23) may lead to a quantitative analysis of the splashing event using the corresponding splashing models.

## V. CONCLUSION

The local limiting crest configuration of progressive solitary interfacial waves confined by two walls parallel to the undisturbed interface is fully determined by three dimensionless numbers: the depth ratio, the density ratio, and the Froude number. In this study, the density ratio is set equal to the small value of 0.001 that corresponds roughly to that of the air-water situation at sea level; the depth ratio is fixed at 80—it is large enough to allow the solitary waves to develop freely. The irrotational flows in both layers are numerically solved using the boundary integral equation method and Newton's iteration method. The solution branch in the amplitude-Froude number plane is obtained, and a highly accurate solution of the limiting configuration with the maximum residual error less than $10^{-11}$ is presented.

Along the solution branch, there exist two different stages. In the first stage, the wave has a non-overhanging crest and the energy increases. Then, a limiting configuration characterized by a droplet sitting on the wave crest with an angle of 120° can be observed. The formation of the droplet-type crest is accompanied by a remarkable drop of both kinetic and potential energies of the water in the crest. This limiting configuration experiences intense shear at the tip, and the shear velocity can be more than five times larger than the wave speed. Furthermore, a simple set of scaling relations is also presented for the fall of the droplet, including the length scale of the droplet and the scales of its impacting velocity and angle. Our results confirm that the droplet-type crest exists even for the low air-water density ratio, without the presence of surface tension. It may be helpful for understanding the common, yet unexplained, phenomena of water waves with a droplet on their crest. The stability analysis is left for future work.

## ACKNOWLEDGMENTS

B. C. Nie would like to acknowledge the support of the National Natural Science Foundation of China (Grant No. 11902024). F. Dias would like to acknowledge the support of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 833125-HIGHWAVE). B. C. Nie and X. Guan would also like to acknowledge the support from the Chinese Scholarship Council.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Bingchuan Nie:** Conceptualization (equal); Writing – original draft (equal). **Xin Guan:** Software (equal); Writing – review & editing (equal). **Jean-Marc Vanden-Broeck:** Writing – review & editing

(equal). **Frédéric Dias:** Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

¹E. Villermaux and Y. Pomeau, “Super free fall,” J. Fluid Mech. 642, 147–157 (2010).

²C. J. Amick, L. E. Fraenkel, and J. F. Toland, “On the Stokes conjecture for the wave of extreme form,” Acta Math. 148, 193–214 (1982).

³L. W. Schwartz, “Computer extension and analytic continuation of Stokes’ expansion for gravity waves,” J. Fluid Mech. 62(3), 553–578 (1974).

⁴J.-M. Vanden-Broeck and L. W. Schwartz, “Numerical computation of steep gravity waves in shallow water,” Phys. Fluids 22(10), 1868–1871 (1979).

⁵J. K. Hunter and J.-M. Vanden-Broeck, “Accurate computations for steep solitary waves,” J. Fluid Mech. 136, 63–71 (1983).

⁶C. W. Lenau, “The solitary wave of maximum amplitude,” J. Fluid Mech. 26(2), 309–320 (1966).

⁷H. Okamoto et al., The Mathematical Theory of Permanent Progressive Water-Waves (World Scientific, 2001), Vol. 20.

⁸J. Witting, “On the highest and other solitary waves,” SIAM J. Appl. Math. 28(3), 700–719 (1975).

⁹D. I. Meiron and P. G. Saffman, “Overhanging interfacial gravity waves of large amplitude,” J. Fluid Mech. 129, 213–218 (1983).

¹⁰J. Y. Holyer, “Large amplitude progressive interfacial waves,” J. Fluid Mech. 93(3), 433–448 (1979).

¹¹P. G. Saffman and H. C. Yuen, “Finite-amplitude interfacial waves in the presence of a current,” J. Fluid Mech. 123, 459–476 (1982).

¹²R. E. L. Turner and J.-M. Vanden-Broeck, “The limiting configuration of interfacial gravity waves,” Phys. Fluids 29(2), 372–375 (1986).

¹³R. H. J. Grimshaw and D. I. Pullin, “Extreme interfacial waves,” Phys. Fluids 29(9), 2802–2807 (1986).

¹⁴D. V. Maklakov and R. R. Sharipov, “Almost limiting configurations of steady interfacial overhanging gravity waves,” J. Fluid Mech. 856, 673–708 (2018).

¹⁵C. J. Amick and R. E. L. Turner, “A global theory of internal solitary waves in two-fluid systems,” Trans. Am. Math. Soc. 298(2), 431–484 (1986).

¹⁶D. I. Pullin and R. H. J. Grimshaw, “Finite-amplitude solitary waves at the interface between two homogeneous fluids,” Phys. Fluids 31, 3550–3559 (1988).

¹⁷T. Kataoka, “The stability of finite-amplitude interfacial solitary waves,” Fluid Dyn. Res. 38(12), 831 (2006).

¹⁸X. Guan et al., “A local model for the limiting configuration of interfacial solitary waves,” J. Fluid Mech. 921, A9 (2021).

¹⁹H. Sha and J.-M. Vanden-Broeck, “Two-layer flows past a semicircular obstruction,” Phys. Fluids A 5, 2661 (1993).

²⁰J. G. Herterich and F. Dias, “Extreme long waves over a varying bathymetry,” J. Fluid Mech. 878, 481–501 (2019).

²¹X. Barthelemy et al., “On a unified breaking onset threshold for gravity waves in deep and intermediate depth water,” J. Fluid Mech. 841, 463–488 (2018).

²²M. S. Longuet-Higgins, “The instabilities of gravity waves of finite amplitude in deep water I. Superharmonics,” Proc. R. Soc. London, A 360(1703), 471–488 (1978).

²³M. Tanaka, “The stability of solitary waves,” Phys. Fluids 29(3), 650–655 (1986).

²⁴P. G. Saffman, “The superharmonic instability of finite-amplitude water waves,” J. Fluid Mech. 159, 169–174 (1985).

²⁵M. Tanaka, “The stability of steep gravity waves,” J. Phys. Soc. Jpn. 52(9), 3047–3055 (1983).

²⁶T. Kataoka, “On the superharmonic instability of surface gravity waves on fluid of finite depth,” J. Fluid Mech. 547, 175–184 (2006).