# Analysis of high-energy drop impact onto deep
liquid pool

Hui Wang1<EMAIL_ADDRESS>Shuo Liu1  Annie-Claude Bayeul-Lainé1  David
Murphy2  Joseph Katz3  Olivier Coutier-Delgosha1,4<EMAIL_ADDRESS>1Univ.
Lille, CNRS, ONERA, Arts et Metiers Institute of Technology, Centrale Lille,
UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de
Fériet, F-59000 Lille, France 2Department of Mechanical Engineering,
University of South Florida, Tampa, FL 33620, USA 3Department of Mechanical
Engineering, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD
21218, USA 4Kevin T. Crofton Department of Aerospace and Ocean Engineering,
Virginia Tech, Blacksburg, VA 24060, USA

###### Abstract

The present work is devoted to the analysis of drop impact onto a deep liquid
pool. It is focused on the effects of high-energy splash regimes, caused by
the impact of large raindrops at high velocities. Such cases are characterized
by short time scales and complex mechanisms, and they have thus received
little attention until now. The BASILISK open-source solver is used to perform
three-dimensional Direct Numerical Simulations (DNS). The capabilities of the
octree adaptive mesh refinement techniques enable to capture the small-scale
features of the flow, while the Volume of Fluid (VOF) approach combined with a
balanced force surface tension calculation is applied to advect the volume
fraction of the liquids and reconstruct the interfaces. The numerical results
compare well with experimental visualizations: both the evolution of crown and
cavity, the emanation of ligaments, the formation of bubble canopy, and the
growth of the downward spiral jet that pierces through the cavity bottom, are
correctly reproduced. Reliable quantitative agreements are also obtained
regarding the time evolution of rim positions, cavity dimensions and droplet
distributions through an observation window. Furthermore, simulation gives
access to various aspects of the internal flows (velocity, vorticity,
pressure, energy budget), which allows to better explain the corresponding
physical phenomena. Details of the early dynamics of bubble ring entrapment
and splashing performance, the formation/collapse of bubble canopy, and the
spreading of drop liquid, are discussed. The statistics of droplet size show
the bimodal distribution in time, corroborating distinct primary mechanisms of
droplet production at different stages.

###### keywords:

##  1 Introduction

The impact of raindrops on a deep liquid pool has been extensively studied
since the initial works of Worthington (1883, 1908). Various behaviours are
obtained, depending primarily on the miscible or immiscible character of the
two liquids (Lhuissier et al., 2013; Castillo-Orozco et al., 2016), but also
on their difference of density (Thomson & Newall, 1886; Manzello & Yang, 2002;
Villermaux, 2022), the depth of the receiving volume (Macklin & Hobbs, 1969;
Wang & Chen, 2000; Fedorchenko & Wang, 2004), and the angle of the impact
(Zhbankova & Kolpakov, 1990; Okawa et al., 2008; Gielen et al., 2017; Liu,
2018), to mention only a few parameters.

In the specific case of identical liquids of density $\rho$, with a
$90^{\circ}$ impact of the drop falling in the air on a deep volume of target
liquid, very different phenomena are still observed, when the liquid viscosity
and surface tension with air, the diameter $d$ of the drop and the speed of
the impact $V$ are varied. Schotland (1960) has first reported the primary
effect of the Weber number $We=\rho V^{2}d/\sigma$ on the physics of the
impact, i.e. the ratio of the drop kinetic energy to the energy required to
deform the target liquid surface.

A few years later, Engel (1966, 1967) has provided a detailed description of
various behaviours obtained from waterdrop impact on water pool when $We$ is
varied in the range $3000\sim 20000$. Based on visualizations of the impact
and using white particles in the target liquid and red ink in the impacting
drop, she has shown that the drop impact creates a cavity in the flat free
surface, and subsequently the rise of the target liquid in a cylindrical
shape, from the edges of the cavity. A bubble-thin cylindrical sheet of liquid
is erected at the upper edge of this crown, which eventually necks in and
closes in a bubble dome in the most energetic test cases. After that, the
cavity shallows and a jet forms at the cavity floor, which flows through the
centre of the crown in case it is open or merges with a downward jet coming
from the top of the bubble dome in case the crown has closed. The author
suggests that the cavity may vibrate at its maximal expansion if all the drop
kinetic energy is not yet transformed into potential energy, and she also
assumes that: (i) the maximum possible bubble height is equal to the cavity
diameter, (ii) the pressure evolution below the cavity floor should explain
the formation of the upward jet, (iii) the liquid of the initial drop is
carried by this jet and will eventually form a secondary bubble at the centre
of the cavity.

In the following decades, this analysis has been progressively completed by
additional experiments with increasing visualization capabilities (Rodriguez &
Mesler, 1985; Pumphrey et al., 1989; Rein, 1996), leading to characteristic
laws of cavity size and growth time (Leng, 2001) and an improved understanding
of the phenomena mentioned hereabove (Rein, 1993). Indeed, these studies were
usually focused on a specific mechanism obtained in a small range of $We$ and
Froude number $Fr=V^{2}/gd$, like for example the crater formation and
collapse (Prosperetti et al., 1989), the break-up of the upward jet for $We$
between 5 and 200 (Manzello & Yang, 2002; Castillo Orozco, 2015), the bubble
entrapment under the floor of the cavity and the subsequent rising thin jet
(Thoroddsen et al., 2003; Deng et al., 2007; Hendrix et al., 2016), the
breakup of the crown rim (Deegan et al., 2008), and the air layer trapped
below the impacting drop when the respective surfaces of liquid are deformed
due to the local pressure changes at the early stage of the process (Tran et
al., 2013).

In the most recent works, attention has been especially focused on these first
stages of the impact, where the entrapped air layer under the drop is
ruptured, resulting in the generation of secondary bubbles (Pumphrey & Elmore,
1990; Thoroddsen et al., 2003, 2012). Immediately after the impact, the
contact line between the drop and the receiving liquid, also called the “neck
region”, moves radially at high speed. It results in a liquid sheet, composed
of water coming from the target liquid, which is ejected horizontally outwards
at the base of the contact surface (Weiss & Yarin, 1999; Thoroddsen, 2002;
Josserand & Zaleski, 2003). This ejecta eventually breaks up into ligaments
and droplets that result in the splashing of very fine sprays (Thoroddsen et
al., 2011; Zhang et al., 2012b). For a sufficiently high Reynolds number, this
phenomenon comes earlier than the sheet-like jet that is directly produced by
the drop impact and will later generate the crown formation, as observed by
Zhang et al. (2012a) using fast X-ray phase contrast imaging. The combination
of these two jets was analysed by Agbaglah et al. (2015) who used a
combination of X-ray imaging and axisymmetric simulations to discuss the
interaction between the two jets.

One important output is that this interaction may be strongly related to
instabilities observed in the neck region. Indeed, the base of the ejecta may
become unstable and generate some alternate vortex shedding in the liquid
below the neck region (Thoraval et al., 2012; Castrejón-Pita et al., 2012). A
strong interplay between this von Kármán vortex street and the
inception/break-up of the ejecta sheet is suggested by several recent studies
(Agbaglah et al., 2015; Li et al., 2018), but the mechanisms are still not
elucidated. Non-axisymmetric behaviours have been observed in the vortex
shedding and the ejecta formation processes from the bottom visualizations in
cases of impact on thin water film (Thoraval et al., 2013; Li et al., 2018),
which shows that these phenomena break the axisymmetry at the fine scale, and
may be related to three-dimensional mechanisms. Local streamwise vortices,
generated at the sharp corners at the basis of the ejecta at an early stage of
the drop impact, are mentioned by these authors to possibly explain these non-
axisymmetric effects. Using ultra-fast imaging, an early azimuthal instability
was also detected by Li et al. (2018), which consists of small azimuthal waves
that grow at the edge of the outer line of the contact, before the inception
of the shedding. This clearly poses a new challenge to numerical studies,
which mostly assume axisymmetry (Thoraval et al., 2012; Agbaglah et al.,
2015). The authors point out the need of three-dimensional simulations with
sufficient resolutions of the finest scale to explore the vortex structure in
the neck region and the interaction between the two corners of the base of the
ejecta sheet, which may be a primary mechanism of this instability.

As for the later stages after impact, phenomena occurring specifically at high
Weber numbers, like the formation of ligaments at the edges of the crown, the
ejection of droplets, the closure of the upper rim, and the subsequent
downward jet, have been also investigated by Engel (1966), Deegan et al.
(2007), Bisighini et al. (2010), Castillo Orozco (2015), and most recently by
Murphy et al. (2015), who provided detailed information about the length of
the ligament according to time, and the population of aerosolized droplets.

Figure 1: Various configurations of splashing reported in the literature
according to the $We$ and $Fr$ numbers reproduced from Murphy et al. (2015)
with the authorization of the authors.

An exhaustive review of these previous works devoted to the impact of a drop
on a deep volume of the same low-viscosity liquid has been conducted by Murphy
et al. (2015), as shown in figure 1. Five different regimes have been
identified by authors, depending on the values of the Froude and Weber
numbers. The two directions of variations of the impact speed $V$ and the drop
diameter $d$ are plotted in figure 1, as well as the specific case of
raindrops of different sizes falling at terminal speed (solid black line
labelled “Raindrop TS”).

As can be seen in figure 1, if the drop size and the impact speed increase
(moving typically from the bottom left to the top right on the chart), the
following successive regimes are obtained: (1) slow drop that coalesces with
the liquid volume and generates a vortex ring that moves downwards as the drop
sinks, (2) formation of a cavity together with a surface wave, (3) the crater
formed by increasing the drop diameter and the surface wave forms a crown rim
from which the secondary droplets and ligaments are ejected, (4) for the
highest speeds and the largest raindrop diameters, a large cavity is obtained,
around which extends vertically a very thin liquid film, whose elevation and
radius keep increasing and it eventually collapses in a very short time. At
the upper part of the film, droplets and ligaments are ejected upwards and
outwards.

This last type of behaviour (black diamonds in figure 1) is the most energetic
and only received very limited attention during the past twenty years (Pan et
al., 2008; Bisighini et al., 2010; Sochan et al., 2018; Lherm et al., 2021).
This high-energy regime is characterized by shorter time scales and increased
complexity of the phenomena involved in the splashing, which makes both
experimental and numerical approaches more challenging, and may explain why it
has received less attention until now. However, it is the case likely to
produce the greatest number of aerosol droplets. The present study is focused
on this highly energetic case of large raindrops falling near terminal speed,
which is representative of raindrops at the surface of the ocean.

This configuration has been studied by Murphy et al. (2015) in detail,
primarily as a reference case for a study focused on the influence of oil
slicks and oil dispersants on the impact. High-speed videos have enabled to
characterise the time evolution of the external shape of the cavity, the crown
and its upper rim, while microscopic holography has provided some statistics
of droplet population ejected at different stages of the process. These
previous experiments are used to validate the numerical strategies conducted
in the present study in section 3, which is focused on the same conditions of
drop diameter and speed at the impact. Further attention is given here to the
early-time bubble ring entrapment and splashing near the impact neck, the
internal mechanisms of crown formation as well as its violent closure, the
liquid jets formed successively upward and downward and the subsequent
entrainment of a large air bubble, and the ejection of droplets at the edge of
the upper rim, in various directions that depend on the stage of the process.
Indeed, it was estimated in the experiments that such a large raindrop
produces at least 2000 micro-droplets. However, the details of the production
of these tiny droplets and statistics on their size, speed, and ejection
direction remain to be found, which motivated the use of numerical simulations
in this paper.

Multiphase flow calculations of drop impacts usually consist in laminar
simulations, as the Reynolds number $Re=\rho Vd/\mu$ based on the drop motion
in the air at the impact velocity is below $10^{4}$, so there is no transition
to turbulence, and all motions resulting from the impact are characterized by
short characteristic times that do not involve any turbulence effect. The
challenge of such computations is thus mostly related to the prediction of the
multiphase structure. Early numerical works (Harlow & Shannon, 1967; Ferreira
et al., 1985) based on Marker and Cell and SOLA-VOF techniques of interface
tracking could provide some first predictions of the main features of drop
impact but could not resolve the small-scale mechanisms. After that, Oguz &
Prosperetti (1990) have included the surface tension effects and Morton et al.
(2000) have solved the full Navier-Stokes equations, opening the way to modern
simulations that are based either on interface tracking (level set, front
tracking) or interface reconstruction, using a transport equation for the
volume fractions of one component (VOF method).

The recent simulations of liquid drops impinging onto the liquid surface are
mostly based on the latter category of models. Both axisymmetric simulations
(Morton et al., 2000; Berberović et al., 2009; Ervik et al., 2014; Ray et al.,
2015; Agbaglah et al., 2015; Deka et al., 2017; Fudge et al., 2021) and three-
dimensional simulations (Rieber & Frohn, 1999; Brambilla & Guardone, 2013;
Cheng & Lou, 2015; Shin et al., 2017) can be found in literature, using mostly
DNS approaches. For three-dimensional calculations, whose primary objective is
to capture the non-axisymmetric mechanisms involved in the splashing, efforts
are focused on dynamic refinement techniques (Nikolopoulos et al., 2007;
Brambilla & Guardone, 2015), in order to resolve the multiple interfaces
resulting from the splashing, while ensuring a reasonable grid size. In the
latter study, the authors show that a minimum cell size of about 50 to 100
$\upmu$m is sufficient to capture most of the features of the splash. However,
that conclusion was drawn for Weber numbers around 250, so it may not directly
apply to the high-energy configuration investigated in the present study.

The paper is organized as follows: the numerical methods and problem statement
are described hereafter in section 2, the detailed validation of the numerical
strategies is conducted by comparing the numerical results with the
experimental data in section 3, the analysis of several mechanisms involved in
the splashing is presented in section 4, and the statistics of airborne
droplets is finally analyzed in section 5.

##  2 Numerical methods and flow configurations

###  2.1 Main features of the solver

The drop impact of the gas-liquid system is considered as incompressible flow,
and it solves a system of the mass balance equation, the momentum balance
equation and the advection of one-fluid formulation, called hereafter the
colour function.

$\nabla\bcdot\boldsymbol{U}=0$ (1) $\rho\left(\frac{\p\boldsymbol{U}}{\p
t}+(\boldsymbol{U}\bcdot\nabla)\boldsymbol{U}\right)=-\nabla
P+\nabla\bcdot(\mu\mathsfbi{D})+\rho\boldsymbol{a}+\sigma\kappa\delta_{s}\boldsymbol{n}$
(2) $\frac{\p C}{\p
t}+\nabla\bcdot(C\boldsymbol{U})=C\nabla\bcdot\boldsymbol{U}$ (3)

Where $\boldsymbol{U}$ is the velocity vector, $\rho$ is the density, $P$ is
the pressure, $\mu$ is the viscosity, $\mathsfbi{D}$ is the deformation tensor
whose components are $D_{ij}={\p u_{i}}/{\p x_{j}}+{\p u_{j}}/{\p x_{i}}$ with
$u_{i}$ and $x_{i}$, $i=1$ to 3 the components of $\boldsymbol{U}$ and
position $\boldsymbol{X}$, respectively. $\boldsymbol{a}$ is the body force
along the impact direction. The last term in equation 2 is the surface tension
force, with $\kappa$ the curvature of the interface, $\sigma$ the surface
tension coefficient which is taken as constant in the present study, and the
$\boldsymbol{n}$ the unit vector normal to the interface. This term is zero
everywhere but at the gas/liquid interface, as controlled by the Dirac
function $\delta_{s}$. The colour function $C$ is transported by equation 3.
For incompressible flow, the right-hand term in equation 3 is zero due to
equation 1.

The BASILISK framework is a flow solver developed by Popinet (2015) and
available for use and development under a free software GPL license (Popinet &
collaborators, 2013–2023). It solves the time-dependent
compressible/incompressible variable-density Euler, Stokes, or Navier-Stokes
equations with second-order space and time accuracy. The Momentum-Conserving
Volume-of-Fluid (MCVOF) approach implemented by Fuster & Popinet (2018) is
employed here to simulate the problem of gas-liquid two-phase flow. The colour
function 3 is solved based on a volume fraction advection scheme proposed by
Weymouth & Yue (2010), which exhibits complete mass conservation and makes it
ideal for highly energetic free-surface flows. The interface then can be
represented using Piecewise Linear Interface Capturing (PLIC) VOF method
(Rudman, 1998). For equation 2, the Crank–Nicholson discretization of the
viscous terms is second-order accurate in time and unconditionally stable,
while the convective terms are computed using the Bell-Colella-Glaz (BCG)
second-order unsplit upwind scheme (Bell et al., 1989), which is stable for
CFL numbers smaller than one. The calculation of the surface tension is one of
the most challenging steps of the process, since no continuous definition of
the interface is available in the VOF approach. Here the balanced-force
surface-tension calculation is used (Francois et al., 2006), which is based on
the Continuum-Surface-Force (CSF) approach originally proposed by Brackbill et
al. (1992). In addition, a second-order accurate calculation of the curvature
is performed, using the Height-Function technique developed by Popinet (2009).
More detailed descriptions of the numerical schemes can be found in Fuster &
Popinet (2018), Pairetti et al. (2020) and Zhang et al. (2020).

Cubic finite volumes organized hierarchically in octree are used for space
discretization. The octree structure has been developed initially for image
processing (Samet, 1990) and later applied to CFD (Khokhlov, 1998) and
multiphase flows (Popinet, 2003). The basic organization is the following, all
details can be found in Popinet (2015): when a cell is refined, it is divided
into 8 cubic cells, whose edges are half the ones of the parent cell. The base
of the tree is the root cell and the last cells with no child are the leaf
cells. The cell level is the number of times it is refined, compared with the
root cell (level 0). To avoid too much complexity in the gradient and flux
calculations, the levels of direct and diagonally neighbouring cells are
constrained and cannot differ by more than one and all cells directly
neighbouring a mixed cell must be at the same level.

BASILISK employs the Adaptive Mesh Refinement (AMR) (van Hooft et al., 2018)
to adaptively refine/coarsen the grids based on the wavelet-estimated
numerical error of the local dynamics, which makes it especially appropriate
for the present application where multiple interfaces and numerous droplets
and bubbles are expected. The resolution will be adapted every time step
according to the estimated discretisation error of the spatially-discretized
fields (volume fraction, velocity, curvature, etc). The mesh will be refined
as long as the wavelet-estimated error exceeds the given threshold, which
eventually leads to a multi-level spatial resolution from the minimum
refinement level $L_{min}$ to the maximum refinement level $L_{max}$ over the
entire domain.

The BASILISK code has been already applied to the study of drop impact onto
liquid film (Josserand & Zaleski, 2003; Reijers et al., 2019; Wu et al., 2021;
Fudge et al., 2021; Sanjay et al., 2022) and the great superiority of parallel
capacity and computational efficiency of the solver were discussed in (Wu et
al., 2021). Furthermore, the capabilities of the BASILISK solver have been
extensively validated on various problems of multi-phase flows here Popinet &
collaborators (2013–2023).

###  2.2 Initial Flow configurations

\begin{overpic}[width=433.62pt]{fig/Fig2.pdf} \put(2.0,55.0){\small({a})}
\put(60.0,63.0){\small({b})} \put(60.0,30.0){\small({c})} \end{overpic} Figure
2: Initial numerical configuration of the three-dimensional simulation. (a)
Overall view of the computational domain and the initial mesh structure at a
plane across the centre of the impacting drop (z = 0). (b) Closeup view of the
initial flows around the impacting drop. (c) Mesh refinement strategy at the
initial stage ($S1$). A higher maximum level of refinement ($L_{max}=14$) is
imposed near the neck region to capture the early-time splashing.

The configuration of drop impact investigated in the present study mimics the
one studied by Murphy et al. (2015), using high-speed video. In the
experiments, the drop falls in a $15.2\times 15.2$ ${\rm cm^{2}}$ tank filled
with seawater to a depth of 8 cm. The measured horizontal and vertical
diameters of the drop just before the impact are $d_{h}=4.3$ mm and
$d_{v}=3.8$ mm, respectively, which results in an effective drop diameter
$d=(d_{v}d_{h}^{2})^{1/3}=4.1$ mm. The ratio between the width of the tank $L$
and the drop diameter $d$ is therefore defined as $L=38d$. The speed before
impact is $U_{0}=7.2$ m/s, which is 81% of the drop terminal speed.

In our numerical simulation, the computational domain is reduced down to a
cube with a side length $L=16d$, which represents about 1/14 of the water tank
volume in the experiments, as shown in figure 2(a). It has been found to be
the best compromise to avoid any effect of the boundary conditions on the
splashing, while decreasing as much as possible the dimensions of the domain.
The free surface is located at the mid-distance between the bottom and the
top, thus the depth of the pool is $H=8d$, to give enough space to the
aerosolized droplets without any interaction with the boundaries. The free
outflow boundary condition is imposed on the top of the domain, while the
default slip boundary condition (symmetry) is applied for the four side walls
and the bottom. A zoomed-in view of the initial flow set-ups in the vicinity
of the drop is depicted in figure 2(b). The initial gap between drop and pool
is $\delta=0.1d$, allowing the observation of air sheet entrainment near the
contact line. The water in the drop and the tank are the same, which is
assumed to have almost no effect on the splashing, as the properties of the
freshwater drop and the target seawater in the experiments are very close. The
density and viscosity ratios between seawater and air are
$\rho_{w}/\rho_{a}=1018$ and $\mu_{w}/\mu_{a}=180$, which leads to a system of
drop-pool impact with dimensionless numbers $Re=28800$, $We=2893$ and
$Fr=1322$.

As illustrated in figure 2(a), the initial mesh configuration around the drop
is generated based on the AMR algorithm using the estimated-discretization
error of volume fraction ($f_{Err}=1e-6$) and velocity ($u_{Err}=1e-4$)
fields, which promises a rounded geometric description of the initial drop and
lowers the RAM requirement of initialization. The mesh is coarsened gradually
down to the given minimum level of refinement ($L_{min}$) away from the drop
interface. The region around the free surface of the pool is refined at
$L_{max}=11$ to avoid any divergence issue. Once the simulation starts, the
mesh will be redistributed adaptively based on AMR using the volume fraction
field with tolerance $f_{Err}=1e-4$ and the velocity field with tolerance
$u_{Err}=1e-1$ as adaption criteria. Additionally, we remove the droplets that
approach the boundaries of the computational domain, since these tiny droplets
have few effects on the evolution of the main impact but are expensive to
track. In real experiments, they are considered to evaporate at one point and
this is not our focus in this paper. Using the initial contact centre as the
reference, any microdroplet whose centroid lies outside of the region of a
semi-sphere with diameter $D_{remove}=15d$ will be removed from the
simulation. The effect of gravity is taken into consideration in this study.

###  2.3 Minimum spatial resolution

Direct Numerical Simulation (DNS) consists in resolving all scales of the
flow, from the smallest relevant ones (ideally here the smallest water
droplets ejected in the air or air bubbles entrained in the water) up to the
large scales of the problem (here the volume of water tank that receives the
initial impacting drop). Ideally, it means that the grid resolution of the
air/water domain should be fine enough to capture all air/water interfaces
created at all steps of the splashing process. We have carried out extensive
tests to explore the effects of the minimum spatial resolution in comparison
with the available experimental data. It was found that the early dynamics of
liquid sheet and air entrapment in the neck region affect significantly the
subsequent phenomena, such as the formation of the crown and its closure at
high $Re$ and $We$ numbers. Our preliminary study has shown that a grid
resolution of at least 1024 cells per drop diameter is necessary to capture
the correct regime of early splashing for the present work (see Appendix A for
details). This will need to apply a maximum level of refinement at
$L_{max}=14$, which corresponds to an equivalent uniform grid of more than 4.3
trillion $[{(2^{14})}^{3}]$ cells. Note that the more the maximum refinement
level is increased, the more the time steps are reduced to comply with the CFL
(Courant–Friedrichs–Lewy) condition, which means that the calculation CPU time
is not proportional to the number of cells. Although the total number of cells
reduces significantly by employing the AMR algorithm, it is still far too
expensive to perform a long-term simulation in full 3D configuration at this
level.

Therefore, based on the physical characteristics during the impact, we divide
the simulation into three consecutive stages, namely the early-time splashing
stage at $t\leqslant 0.2$ ms ($S1$), the crown formation stage at $0.2<t<4$ ms
($S2$) and the bubble canopy stage at $t\geqslant 4$ ms ($S3$). At stage $S1$,
very thin liquid sheet emerges on the neck region and interacts strongly with
the surface of drop and pool, producing numerous very fine droplets and
bubbles, thus the primary objective is to capture the flow dynamics near the
contact region in finer scales. Figure 2(c) shows the mesh refinement strategy
employed at the initial stage. A higher maximum level $L_{max}=14$ (3.9
$\upmu$m) is used in the vicinity of the free surface to solve the dynamics of
the neck region in smaller scales, while $L_{max}=13$ (7.8 $\upmu$m) is used
for the rest of the domain to capture the general dynamics. At stage $S2$, a
coherent liquid sheet has been developed above the pool and it grows
subsequently into the thin-walled crown. During this rapid expansion stage, a
great deal of secondary droplets is sprayed incessantly from the top of the
crown (Murphy et al., 2015). The extra layer of refinement with $L_{max}=14$
is thus removed at this stage and $L_{max}=13$ is used for the entire
computational domain to capture the statistics of droplets and bubbles. At
stage $S3$, the droplets shed from the top of the crown are much less and
generally larger, comparing with the ones from previous stages. Therefore, we
restart the simulation with $L_{max}=12$ (15.6 $\upmu$m) over the whole
computational domain, which allows capturing the main physical dynamics of
crown and cavity till the end of the simulation (t = 48 ms). This mesh
refinement strategy ensures that the simulation is doable in a full three-
dimensional configuration and can be accomplished in a “reasonable” time,
considering the available computational resources for the moment.

All the numerical results presented in the main text of this paper were
performed on 1024 cores for 33.5 days, which consumes more than
$8.21\times{10}^{5}$ CPU-hours in total, using the computational resources on
Advanced Research Computing (ARC) at Virginia Tech.

##  3 Comparisons with the experiments

###  3.1 Morphology

\begin{overpic}[scale={1.2}]{fig/Fig3.pdf} \put(-4.0,62.4){\small({a})}
\put(-4.0,22.2){\small({b})} \end{overpic} Figure 3: Qualitative comparison
between experiments and simulations. (a) Overall dynamics of the air-water
interface during 48 ms after impact. From left to right, the experiment shows
-1, 1, 3, 7, 12, 18, 41 and 52 ms after impact, and the simulation shows -1,
1, 3, 7, 12, 18, 37 and 48 ms after impact. The red stars indicate the tracked
positions of the upper rim of the crown. The scale bar is 10 mm long. (b)
Closeup view of the early-time splashing during 450 $\upmu$s after impact.
From left to right, the experiment shows 49, 148, 246, 345 and 443 $\upmu$s
after first contact, and the simulation shows 50, 150, 250, 350 and 450
$\upmu$s after first contact. The scale bar is 1 mm long. The qualitative
comparisons show that the simulation successfully reproduced all the
distinctive features observed in the experiments.

Figure 3(a) shows the global evolution of the air-water interfaces generated
by high-energy drop-pool impact during the 48 ms after contact. The side-views
of the numerical results (bottom) are compared with the experimental high-
speed images (top) in time. The time point of contact is defined as $t=0$ ms.
At $t=-1$ ms, the drop is initialised above the surface of the pool. Once the
simulation starts, the drop will fall downwards to hit the liquid surface
driven by the combination of the initial impact speed and the gravitational
force. Directly after impact ($t=1$ ms), a cylindrical wave around a flat-
bottomed disk-like cavity is produced due to the drop violent penetration and
substantial liquid ligaments emanate almost horizontally from the thickened
rim of the crown, spraying a large number of micro-droplets into the air. At
this stage, most of the drop liquid is concentrated on the bottom of the
cavity and the submergence of the drop/pool interface can be approximated as
half of the initial drop impact velocity (Fedorchenko & Wang, 2004). In the
following few milliseconds, the drop keeps expanding radially and eventually
spreads into a thin layer of liquid that is distributed along the interior
surface of the cavity, stretching the cavity into a typical hemispherical
shape that has been widely reported in literature (Engel, 1966; Berberović et
al., 2009; Bisighini et al., 2010; Ray et al., 2015), as seen at $t=3$ ms. By
$t=7$ ms, the upper rim of the crown has started to proceed inwards and the
orientation of the ligaments transitioned from almost horizontal to vertical,
which eventually leads to the upcoming closure event on the upper part. At the
time when the crown necks in ($t\approx 12$ ms), flows from the sidewalls of
the cylindrical wave meet along the impact axis and a large volume of air is
encapsulated, generating a central liquid jet that moves spirally from the
merging point ($t=18$ ms). The downward-moving jet keeps growing and
eventually pierces the cavity bottom as shown at $t=37$ ms, disturbing the
retraction of the cavity bottom. Finally, the rebound of the cavity produces a
broad upward jet that merges with the previous central column of fluid,
leaving several air bubbles inside the liquid tank, evidenced at $t=48$ ms.

Furthermore, a close-up view of the early-time splashing near the contact
region during the first 450 $\upmu$s after impact is provided and compared
with the high-speed experimental holograms in figure 3(b). An immediately
ruptured “liquid sheet” appears right after contact ($t=50$ $\upmu$s), known
as the “prompt splash” (Deegan et al., 2007). In this process, a great number
of very fine droplets are scattered in the air, which has been considered as
one of the primary sources of marine aerosols and may raise potential health
issues to the public under certain circumstances (Murphy et al., 2015). More
liquids are subsequently pushed out from the pool to form the liquid wall of
the crown. The morphological behaviours of the early-time splashing captured
by simulation are quite consistent with the experimental observations.

These qualitative comparisons show that our simulation reproduced all the
distinctive features observed in the experiments. The correct prediction of
the early-time splashing, the transition of the ligaments’ orientation, and
the exact time when the upper rim of the crown necks in and the central spiral
jet pierces the cavity bottom are especially convincing. A general good
agreement to the experiments is obtained.

###  3.2 Kinematics of crown and cavity

\begin{overpic}[width=433.62pt]{fig/Fig4.pdf} \put(0.0,97.5){\small({a})}
\put(35.0,97.2){\tiny{Experiment}} \put(35.0,95.7){\tiny{Simulation}}
\put(21.7,67.5){\scriptsize{$t$ $\rm{(ms)}$}}
\put(1.0,82.0){\scriptsize\rotatebox{90.0}{{$R_{r}$ $\rm{(mm)}$}}}
\put(0.0,64.5){\small({b})} \put(35.0,64.2){\tiny{Experiment}}
\put(35.0,62.5){\tiny{Simulation}} \put(21.7,34.5){\scriptsize{$t$
$\rm{(ms)}$}} \put(1.0,50.0){\scriptsize\rotatebox{90.0}{{$R_{h}$
$\rm{(mm)}$}}} \put(11.0,31.0){\small({c})} \put(23.0,30.5){\tiny{Experiment}}
\put(23.0,29.0){\tiny{Simulation}} \put(20.0,1.0){\scriptsize{$R_{r}$
$\rm{(mm)}$}} \put(12.0,15.0){\scriptsize\rotatebox{90.0}{{$R_{h}$
$\rm{(mm)}$}}} \put(43.8,97.5){\small({d})} \put(79.0,97.2){\tiny{Experiment}}
\put(79.0,95.7){\tiny{Simulation}} \put(65.5,67.5){\scriptsize{$t$
$\rm{(ms)}$}} \put(45.0,82.0){\scriptsize\rotatebox{90.0}{{$C_{r}$
$\rm{(mm)}$}}} \put(43.8,64.5){\small({e})} \put(63.0,64.2){\tiny{Experiment}}
\put(63.0,62.5){\tiny{Simulation}} \put(63.0,61.0){\tiny{Analytical prediction
of \cite[citet]{\@@bibref{Authors
Phrase1YearPhrase2}{bisighini2010crater}{\@@citephrase{(}}{\@@citephrase{)}}}}}
\put(65.5,34.5){\scriptsize{$t$ $\rm{(ms)}$}}
\put(45.0,50.0){\scriptsize\rotatebox{90.0}{{$C_{d}$ $\rm{(mm)}$}}}
\put(43.8,31.0){\small({f})} \put(80.0,30.5){\tiny{Experiment}}
\put(80.0,29.0){\tiny{Simulation}} \put(66.5,1.0){\scriptsize{$t$
$\rm{(ms)}$}} \put(44.5,16.0){\scriptsize\rotatebox{90.0}{{$C_{v}$
$\rm{(mm)}$}}} \end{overpic} Figure 4: Analysis of the quantitative data
measured with respect to the initial impact centre. (a) Evolution of the crown
radius. (b) Evolution of the crown height. (c) Trajectory of the upper rim of
the crown. (d) Evolution of the cavity radius. (e) Evolution of the cavity
height. (f) Evolution of the cavity volume. The black dashed line in (e) shows
the theoretical prediction of the penetration depth using the proposed model
by Bisighini et al. (2010).

The kinematic behaviours of the crown and the cavity are manually tracked from
experiments and simulations, and compared quantitatively in figure 4. The
impact centre where the drop firstly contacts the pool is used as the
reference point. The time evolution of the upper rim of the crown (marked by
red stars in figure 3), its radial distance (figure 4a) and height (figure
4b), are measured during 24 ms after impact. The height is defined as the
distance between the initial quiescent free surface of the pool and the
position where ligaments are forming. The trajectory of the upper rim of the
crown is then plotted in figure 4(c). For the first few milliseconds, the
crown expands rapidly along the horizontal radial direction and reaches its
maximum radius very soon ($t\approx 3$ ms). After the maximum horizontal
position, the rim starts to proceed almost immediately inwards, primarily
under the effect of surface tension. Meanwhile, the leading edge of the crown
rises continuously in the vertical direction and approaches its maximum height
right before it necks in ($t\approx 12$ ms). After the closure of the upper
part, this point vibrates slightly near the impact axis along with the
shrinking toroidal large air bubble.

The evolution of the submerged cavity has been thoroughly discussed in
experiments, simulations and theories, and substantial attention has been
particularly given to the estimation of the geometric dimensions of the cavity
in previous studies (Engel, 1967; Prosperetti & Oguz, 1993; Leng, 2001;
Berberović et al., 2009; Bisighini et al., 2010; Jain et al., 2019). Figure
4(d) demonstrates the temporal variation of the horizontal radius at the
intersection edge between the cavity and the initial free surface. The cavity
grows rapidly on the horizontal plane at the beginning due to the initial
violent impingement and the outward expanding speed decreases gradually. After
the closure of the crown, the increase of the cavity radius slows down and is
nearly linear. Starting from $t\approx 3$ ms, the horizontal radius at the
bottom of the crown becomes wider than at the top, which may presumably
provide an overall momentum towards the central axis, thus propelling more
liquid from the receiving liquid to the sidewalls of the crater. Figure 4(e)
shows the time evolution of the maximum depth of the cavity. The cavity keeps
expanding in depth and reaches its maximum position 24 ms after impact, which
takes almost double the time of the maximum crown position ($t\approx 12$ ms).
The volume of the cavity is calculated as half of an ellipsoid as used in
Murphy et al. (2015) in figure 4(f).

Figure 4 confirms this conclusion: both the trajectory of the edges of the
upper rim of the crown and the geometric dimensions of the cavity are found in
very good agreement with experimental measurements. A very reliable
quantitative agreement is obtained between simulations and experiments.

###  3.3 Droplets

\begin{overpic}[width=433.62pt]{fig/Fig5.eps} \put(1.0,81.0){\small({a})}
\put(22.0,43.0){\scriptsize{$X$ $\rm{(mm)}$}}
\put(2.0,62.0){\scriptsize\rotatebox{90.0}{{$Y$ $\rm{(mm)}$}}}
\put(44.5,81.0){\small({b})} \put(71.5,43.0){\scriptsize{$X$ $\rm{(mm)}$}}
\put(51.0,62.0){\scriptsize\rotatebox{90.0}{{$Y$ $\rm{(mm)}$}}}
\put(1.0,39.0){\small({c})} \put(33.5,38.3){\tiny{Simulation}}
\put(33.5,36.5){\tiny{Experiment}} \put(22.0,1.0){\scriptsize{$Y$
$\rm{(mm)}$}}
\put(2.0,20.0){\scriptsize\rotatebox{90.0}{{$\overline{N}(Y,t)$}}}
\put(44.5,39.0){\small({d})} \put(89.5,38.3){\tiny{Simulation}}
\put(89.5,36.5){\tiny{Experiment}} \put(72.0,1.0){\scriptsize{$S_{d}$
$\rm{(\upmu m)}$}}
\put(46.0,20.0){\scriptsize\rotatebox{90.0}{{$\overline{N}(S_{d},t)$}}}
\end{overpic} Figure 5: Comparisons of the droplet statistics between
numerical and experimental data captured in a specific field of view. (a)
Overall schematic view of the relative position of the observation window. (b)
Closeup view of the secondary droplets in the observation window. (c) Vertical
distribution of secondary droplets. (d) Size distribution of secondary
droplets. The numerical droplet statistics presented in (c) and (d) are time-
averaged statistics using 9 time slices over the time period $3\sim 4$ ms. The
experimental data are ensemble-averaged using more than 25 replications as
originally presented in figure 17 and 18 in Murphy et al. (2015).

The production of secondary droplets and their distribution induced by the
process of normal (perpendicular) (Okawa et al., 2006; Guildenbecher et al.,
2014; Li et al., 2019; Wu et al., 2020) and oblique (Okawa et al., 2008; Liu,
2018) drop impact on a liquid surface can be found in literature, focusing
mostly on the fairly “smooth” splashing under the relatively lower range of
$Re$ and $We$ numbers. The rather violent splashing behaviours observed under
high-energy impacting conditions suggest that the creation of these tiny
droplets might be associated with more complicated “irregular” interfacial
deformation and breakups (Thoraval et al., 2012; Murphy et al., 2015).
Nevertheless, the understanding of the governing mechanisms and their
population remain insufficient, and the available statistical data on droplet
production is very limited in the literature.

Therefore, for further validating the present numerical strategies, the
statistics of secondary droplets are extracted from simulations and compared
with experimental data. The specific location of the observation window is
shown in figure 5(a), where a $10\times 10$ $\rm{mm^{2}}$ square field of view
is placed 13 mm above the free surface of the pool. In the experiments, the
holographic frames for droplet analysis were selected at the time when the
first upward-rising droplet exits the top of the observation window ($t\approx
3\sim 4$ ms), and the droplet statistics are then ensemble-averaged over all
replicates (more than 25), as explained in Murphy et al. (2015). For the
simulation here, the time-averaged droplet statistics are obtained using 9
time slices over the time period $3\sim 4$ ms. Figure 5(c) shows the droplet
distribution in the vertical direction. Similar to the experimental
observations, most of the droplets are still concentrated at a lower position,
which can be clearly seen in figure 5(b). Figure 5(d) plots the time-averaged
droplet diameter distribution. Comparing with the bimodal distribution in
experiments, an inapparent second peak can be still found from the numerical
results under relatively wider bins (doubling the size of bins used in Murphy
et al. (2015)). Two primary plateaux centred around 50 $\upmu$m and 225
$\upmu$m become more pronounced, which is in good agreement with the
experimental data. The successful reproduction of droplet statistics in a
specific observation window is highly encouraging, as it gives confidence for
further comprehensive and in-depth analysis over the broader spatial domain
and distinct temporal stages in the following sections.

##  4 Overall dynamics of splashing

###  4.1 Early-time dynamics

The very early dynamics that occur shortly after contact ($<100$ $\upmu$s) in
the process of drop impact onto liquid surface has been widely discussed
during the past thirty years (Weiss & Yarin, 1999; Deegan et al., 2007; Liang
& Mudawar, 2016; Marcotte et al., 2019). Conventionally, for low impact
velocities, an air disk is entrapped under the centre of the impacting drop by
lubrication pressure and it later breaks up into chains of micro-bubbles
(Thoroddsen et al., 2012; Tran et al., 2013) or contracts into several bubbles
along the central line (Thoroddsen et al., 2003; Jian et al., 2020) due to
instability. As the contact line expands outwards radially, sheet-like liquid
ejecta is sent out nearly axisymmetrically from the outer edge of the neck and
possibly emits rings of tiny droplets from its rim for sufficiently large
Reynolds numbers (Weiss & Yarin, 1999; Thoroddsen, 2002; Josserand & Zaleski,
2003; Howison et al., 2005; Deegan et al., 2007; Marcotte et al., 2019).
However, with increased impact velocity ($Re>7000$), azimuthal undulations are
found to grow at the base of the ejecta and the entrapment of bubble rings
takes place near the neck region (Thoraval et al., 2013; Li et al., 2018),
which thereby breaks the axisymmetry of the motions in fine scales. For even
higher impact velocities, irregularly distributed splashing and incoherent
liquid sheets are experimentally observed by Thoroddsen (2002). Axisymmetric
simulations of Thoraval et al. (2012) have suggested that the base of the
ejecta may become highly unstable under high-energy configurations, which in
turn propels the ejecta base to swing to collide with drop and pool, trapping
alternate air volume on both sides of the ejecta. The existence of the von
Kármán-type vortical structure was experimentally confirmed later by
Castrejón-Pita et al. (2012) using Shadowgraph imaging and laser-sheet
visualizations. Nevertheless, the study of such complex flow structures is
still very challenging for both experimental and numerical approaches as it
mainly occurs in microscopic length within the time scale of several
microseconds after contact.

\begin{overpic}[width=433.62pt]{fig/Fig6.eps} \put(0.0,35.0){\small({a})}
\put(21.5,1.0){\scriptsize{$(t-t_{0})U_{0}/R_{d}$}}
\put(0.5,18.0){\scriptsize\rotatebox{90.0}{{$R_{n}/R_{d}$}}}
\put(49.5,35.0){\small({b})} \put(75.2,1.0){\scriptsize{$tU_{0}/d$}}
\put(51.0,19.5){\scriptsize\rotatebox{90.0}{{$\theta(^{\circ})$}}}
\end{overpic} Figure 6: Early-time dynamic behaviours of the neck region. (a)
Evolution of the neck radial position $R_{n}$. The solid line shows the
theoretical estimate using the form $R_{n}\sim C\sqrt{3(t-t_{0})U_{0}/R_{d}}$,
where $C=1.22$ and $t_{0}=-2.22$ are obtained by fitting the numerical
measurements. (b) Evolution of the ejecta angle $\theta$ measured from the
vertical central slices as shown in figure 8 (two sides), using the definition
sketch proposed by Thoraval et al. (2012). The sharp decrease of the ejecta
angle indicates the “bumping” event, where the contact point of the neck
suddenly changes due to the reconnection between the ejecta and the drop/pool.

Upon contact, most of the liquid remains unperturbed while the connection of
two liquid masses propagates outwards instantaneously. Theoretically, the
spreading law of $R_{n}\sim\sqrt{2\tau}$ can be derived from the truncated
sphere approximation based only on the geometric considerations (Rioboo et
al., 2002; Josserand & Zaleski, 2003), but it violates the continuity
equation. Combining Wagner’s theory (Wagner, 1932), the radial motion of the
neck for drop impact on a solid surface can be described using the form
$R_{n}\sim C\sqrt{3\tau}$ (Riboux & Gordillo, 2014; Philippi et al., 2016; Li
et al., 2018). The dimensionless time is defined as
$\tau=(t-t_{0})U_{0}/R_{d}$, where $R_{d}=d/2$. Figure 6(a) confirms the good
match between simulated neck motion with analytical prediction. The
coefficient $C=1.22$ and $t_{0}=-2.22$ are obtained by fitting the numerical
measurements. At high impacting $Re$, the initial spreading speed of the neck
may reach as high as $\sim 17$ times the impact velocity (Li et al., 2018),
which basically requires the time resolutions of the order of $10^{-7}\sim
10^{-6}$ s to observe the instantaneous motions along the contact line,
possibly explaining the limited knowledge in this specific area.

Figure 7: Early-time dynamics near the contact region induced by high-speed
drop impact onto deep liquid pool observed from the bottom view. The first
three frames are shown 4, 6 and 8 $\upmu$s after the first contact, where the
“nearly axisymmetric” bubble rings are entrapped from the neck of the
connection. The black arrow points at the central air disk. The red arrows
indicate the formation of a new bubble ring. The last four images show a
smaller magnification 10, 15, 32 and 50 $\upmu$s after impact. The outer edge
is the downward-moving drop, the inner edge is the contact line of the neck
and the central irregular disc is the entrapped air sheet. Azimuthal
instabilities and liquid ejecta are developed along the neck. The outer line
of the neck has not reached the size of the impacting drop here. The scale bar
is 500 $\upmu$m long.

Figure 7 shows the simulated early-time dynamics in the neck region between
drop and pool from the bottom view. Here it is focused on the stage before the
outer edge of the neck overtakes the outline of the drop. The first frame and
the last frame reach $25\%$ and $89\%$ of the drop size respectively. The air-
water interface is coloured by the volume fraction of the passive tracer (see
section 4.5) and the opacity of the interface is set to 0.5, which enables to
visualize the interfacial dynamics for both the inside and outside regions of
the neck.

Immediately after the contact, a central disc of air (black arrow) is
entrapped and an outer liquid edge connecting the drop and pool is formed. The
contact line then expands rapidly in the radial direction, entrapping bubble
rings along the drop/target boundary. In the second panel of figure 7, it can
be seen that a new round of air cylinder (red arrow) is entrapped at the neck
and is pinched off shortly into the bulk within the time scale of $<0.4$
$\upmu$s. By $t=8$ $\upmu$s, up to 10 bubble rings are present within the
target volume and most of them are entrapped axisymmetrically as concentric
circles. For the next $\sim 20$ $\upmu$s, these intact air rings will be
stretched longitudinally while rotating, and eventually break up into a
necklace of micro-bubbles due to surface-tension Rayleigh instability
(Chandrasekhar, 2013) as shown at $t=32$ $\upmu$s and $t=50$ $\upmu$s, which
looks very similar to the experimental observations of figure 7 in Thoroddsen
et al. (2003).

\begin{overpic}[scale={0.84}]{fig/Fig8.pdf} \put(-3.5,37.0){\small({a})}
\put(-3.5,29.5){\small({b})} \put(-3.5,19.5){\small({c})} \end{overpic} Figure
8: Flow field and vorticity structure in the vicinity of the neck region on
the vertical slice at $z=0$ (see dashed line at $t=10$ $\upmu$s in figure 7).
The red and blue colours represent counterclockwise and clockwise rotation
respectively. (a) Entrapment of air disk and bubble rings 4 $\upmu$s after the
first contact. The scale bar is 50 $\upmu$m long. (b) Formation of liquid
ejecta from the neck at $t=12$ $\upmu$s. Secondary droplets are emitted from
its tips. Vortex shedding of the alternate sign from the base of the ejecta
generates a Von Kármán-type structure along the drop/pool boundary, with only
occasional air bubbles/bubble arcs entrapment from the neck. The scale bar is
100 $\upmu$m long. (c) Collisions between the ejecta and the downward-moving
drop, so-called the “bumping event”, which leads to the entrapment of air
rings at $t=73$ $\upmu$s. The scale bar is 500 $\upmu$m long.

Note that the neck of connection between drop and pool remains rather smooth
at this time ($t<10$ $\upmu$s) and no valid liquid ejecta is observed around
it, implying that the fundamental mechanism of bubble ring entrapment here
differs from the jet-induced air encapsulation predicted by Weiss & Yarin
(1999). Figure 8 shows the air-water interface overlaid by the vorticity field
near the neck region on the vertical intersection across the drop centre. An
air void is rolled up along the entire circumference of the neck and later
entrapped axisymmetrically inside the liquid phase. This is reminiscent of the
entrapment of toroidal bubbles numerically captured by Oguz & Prosperetti
(1989) for two drops collision, where the basic question “Does the contact
line between two approaching surfaces move outwards fast enough to prevent
further contact after the initial one? ” was discussed. At the moment of
contact, a very small liquid bridge of radius $R_{n}$ with high curvature
“meniscus” connects two liquid masses and a thin outer air sheet retracts
outwards rapidly driven by surface tension and local pressure gradient. A
self-similar repeated reconnection of this air gap was predicted at the very
early time of contact, thus enclosing a number of tiny toroidal bubbles in the
liquid phase. Focusing on the viscous regime of $R_{e}\ll 1$, $R_{e}=\sigma
R_{n}/(\rho\nu^{2})$, the analytical analysis of Eggers et al. (1999) draws
the scaling laws for the radius of the entrapped toroidal bubble $r_{b}\propto
R_{n}^{3/2}$ and the width of the thin air gap connecting the bubble
$r_{g}\propto R_{n}^{2}$. Further investigation of Duchemin et al. (2003) for
inviscid drop coalescence explained that the occurrence of each pinch-off
event ($i_{th}$) depends only on the local dynamics of air gap width
$r_{g}^{i}$, where $r_{g}^{i}=(R_{n}^{i})^{2}$. The distance between the
initial tip of the meniscus and the reconnection point as well as the time
interval can be estimated as $r_{c}^{i}=10(R_{n}^{i})^{2}$ and
$t_{c}^{i}=7.6(R_{n}^{i})^{3}$ respectively. The author emphasised that the
reconnection ceases when $R_{n}>0.05$, thus $r_{c}^{i}\leqslant 0.025$ and
$t_{c}^{i}\leqslant 0.00095$, where the time and space coordinates are scaled
by $\sqrt{\rho{R_{d}}^{3}/\sigma}$ and $R_{d}$. The evolution of the bubble
rings was however not able to be predicted by the analytical model because of
their highly non-circular shape and three-dimensional rotations and
stretching. Based on this theoretical model, the distances and time intervals
of the neighbouring bubble rings for the present case should be $\leqslant
51.25$ $\upmu$m and $\leqslant 10.41$ $\upmu$s. Measurements from our
simulation show that the maximum distance and time interval between bubble
rings are around $\sim 31$ $\upmu$m and $\sim 0.5$ $\upmu$s, in the same order
as the analytical estimate. It should be noted that the existence of the
entrapped central air disk as well as the variation of the drop bottom
curvature at the moment of impact may alter the local dynamics, but the
phenomena should be qualitatively similar.

Figure 9: Irregular azimuthal undulations on the neck region between drop and
pool 10 $\upmu$s after the first contact. The central irregular plate is the
entrapped air disk. The white circle indicates the early-time breakups of the
ejecta fingers. The scale bar is 100 $\upmu$m long.

Starting from $t\approx 10$ $\upmu$s, an azimuthal instability is observed
along the neck. Here the diameter of the outer contact line reaches only
$38\%$ of the drop size. Figure 9 shows the closeup view of the azimuthal
undulations at the neck of connection between drop and pool at $t=10$
$\upmu$s, where the entire periphery of the neck is visible. Irregular
undulations are captured on both sides of the neck, namely the outer liquid
ejecta and the inner air sheet. Liquid fingers are initiated arbitrarily from
the outer side and some of them break up immediately on their tips (white
dashed circle), producing the very first generation of secondary droplets (see
also figure 8b). At some locations where fingers are not found (or long
wavelength), smooth ejecta and air sheets are present. We have measured that
the characteristic wavelength and amplitude of the initial fingers are around
$\sim 30$ $\upmu$m and $\sim 21$ $\upmu$m, which are generally longer than the
thickness of the ejecta base ($\sim 7$ $\upmu$m). As the neck spreads
radially, alternate vortices are shed from the base of the ejecta, thus
forming a similar structure of von Kármán-type street along the drop/target
contact as shown in 8(b). Meanwhile, the oscillations of the base of the
ejecta pull the local air sheet on the corner of it, entrapping occasionally
some isolated bubbles/bubble arcs on both/alternate sides of the base (green
arrow in figure 7).

After the initial small amplitude oscillations combined with the “complete”
disintegration of the ejecta jets, a coherent ejecta sheet eventually emerges
(see $t=32$ $\upmu$s in figure 7). This cylindrical ejecta sheet then rises
radially along the contact line until it impacts the drop surface, entrapping
a larger bubble ring along the neck as evidenced at $t=50$ $\upmu$s in figure
7 and figure 8(c). Figure 6(b) plots the early-time evolution of the ejecta
angle $\theta$ measured based on the definition sketch proposed in Thoraval et
al. (2012). It can be found that $\theta$ increases almost linearly at the
initial time and suddenly decreases at the moment of “bumping”, which has been
also observed by Thoroddsen et al. (2011) and Thoraval et al. (2012). Once the
new neck of connection is established, the fast-moving rim will stretch the
ejecta sheet and tear it immediately into multiple liquid “tori”, which later
break up and produce a large number of similar-sized microdroplets, as
illustrated in figure 8(c) and discussed later in section 5.

###  4.2 Formation of bubble canopy

\begin{overpic}[width=433.62pt]{fig/Fig10.pdf}
\put(21.4,3.4){\scriptsize{$U/U_{0}$}} \put(75.2,3.4){\scriptsize{$P$}}
\end{overpic} Figure 10: Internal flows of the high-speed drop impact
overlapped by velocity and pressure fields. From top to bottom and left to
right, the corresponding times are 1, 3, 11, 16, 24, 32, 38 and 48 ms after
impact. The velocity magnitudes are scaled by the drop impact speed $U_{0}$.
The scale bar is 10 mm long.

Figure 10 shows the internal flows of the liquid phase at different stages,
overlapped by velocity vector and pressure fields on the cross-section. The
grid-based arrow is oriented by the velocity field and its value is
represented by colour. Within the first millisecond, the drop has not yet
fully sprawled out and most of the momentum is still concentrated in the
impacting drop, thus generating a broad high-pressure area along the drop/pool
interface. Meanwhile, an outward-expanding liquid sheet arises along the
contact line as the result of the violent extrusion. As the drop moves
downwards, more liquid is therefore pushed away from the pool and transported
to the uprising crown. At high impact velocities, thin liquid ligaments are
generated along the rim of the crown and break up on its tips by instability,
producing moderate and large-scale secondary droplets (see section 5). As
discussed above in section 3.2, the crown reaches a maximum horizontal radial
position soon after impact ($t\approx 3$ ms) and then its rim thickens and
bends towards the impact axis under the effect of surface tension, while it
simultaneously rises in the vertical direction. Moreover, the outline of the
cavity expands rapidly on the surface of the pool and overtakes the crown
expansion at around $t\approx 4$ ms, which also facilitates generating an
inward-directed momentum on the crown top to enclose the upper part. During
the period of crown expansion, the maximum velocity is reached on the top
where the liquid film is thinnest. Right before the closure ($t\approx 12$
ms), it is visible that the velocity on the crown rim points almost
horizontally inwards. Flows from all directions on the liquid wall meet and
interact at the instant of closing, generating a sharp pressure rise around
the point of closure, which later protrudes upward and downward moving jets
evidenced at $t=16$ ms. Meanwhile, ligaments above the dome tangle and merge,
possibly shedding several large-scale droplets on the top of the dome, as
shown at $t=16$ ms.

It can be found in the literature that large bubble entrapment owing to drop
impact onto liquid pool occurs mainly in two types of mechanisms. Firstly,
under low impact energy, a vortex-induced roll jet may be formed and grows
into a thick liquid tongue, which later collapses near the pool surface to
entrap the air bubble (Pumphrey & Elmore, 1990). The shape of the drop at the
time of impact is crucial for such mechanism and it occurs almost exclusively
for prolate-shaped drops (Zou et al., 2012; Wang et al., 2013; Thoraval et
al., 2016; Deka et al., 2017). Secondly, with sufficiently high impact energy,
the thin-walled liquid crown rises up higher above the pool due to the violent
collision and its rim bends towards the impact axis while rising vertically,
enveloping the air bubble above the target pool. The formation of such complex
flow structures for high-velocity drop-pool collision has been observed from
time to time by experiments during the past century (Worthington, 1908; Engel,
1966; Bisighini et al., 2010; Murphy et al., 2015; Lherm et al., 2021), but is
probably firstly solved in 3D by high-resolution numerical simulations in the
present work. Interestingly, although driven by different mechanisms,
similarities are still observed between them, such as the generation of the
central jets, the reconnection of the downward-moving jet as well as the
bursting of the final “floating bubble”. Comparing with the low-energy
counterparts, the rather intense interfacial deformation at high impact
velocities surely introduces some new phenomena, which could accordingly
influence the subsequent physical process like the generation of underwater
noise induced by air bubble entrainment (Prosperetti et al., 1989; Prosperetti
& Oguz, 1993) and the additional source of airborne droplets caused by liquid
film rupture (Resch et al., 1986; Afeti & Resch, 1990; Resch & Afeti, 1991).
Further investigation and analysis could be initiated in the future to compare
and contrast these seemingly similar dynamic patterns.

\begin{overpic}[scale={1.0}]{fig/Fig11.pdf}
\put(94.0,42.0){\scriptsize{$U/U_{0}$}} \put(-3.0,48.0){\small({a})}
\put(27.2,48.0){\small({b})} \end{overpic} Figure 11: Formation of the central
spiral jet inside the bubble canopy. (a) Dynamics of the liquid jet at $t=24$
ms. The scale bar is 4 mm long. (b) Jet motions observed from the bottom view,
showing 18, 23, 24 and 38 ms after impact. The scale bar is 1 mm long. The
interface is contoured by velocity field.

Figure 11 shows the formation of the central jet from the top of the dome. The
jet is formed side-biased as the flows on the upper part of the crown are not
perfectly axisymmetric and they do not arrive at the merging point
simultaneously, which is also consistent with one of the few available
experimental observations in literature (Bisighini et al., 2010; Murphy et
al., 2015; Lherm et al., 2021). High-speed jet ejections out of a liquid
interface are commonly observed in many other physical processes, such as
bubble bursting (Boulton-Stone & Blake, 1993; Thoroddsen et al., 2009; Berny
et al., 2022), Faraday waves (Hogrefe et al., 1998; Zeff et al., 2000), and
cavity collapse induced by the process of solid/liquid object impact onto
fluid target (Worthington & Cole, 1897, 1900; Ray et al., 2015; Gekle &
Gordillo, 2010; Jamali et al., 2020; Kim et al., 2021). In general, all these
jets are ejected as a consequence of a very large axial pressure gradient
created at the jet base, which therefore can be further classified based on
the way the large overpressure is created and the length scale at which
pressure variations take place (Gekle & Gordillo, 2010). Surface tension plays
an essential role in the jet formation process.

###  4.3 Cavity contraction

Now we focus on the contraction of the cavity. For drop impact at low and
moderate velocities, the liquid crown collapses and generates capillary waves
that move radially towards the cavity bottom along the interior of the crater
after it reaches the maximum expansion. The capillary waves then meet
concentrically at the bottom of the cavity, forming a classic upward
Worthington jet that may break up on its tips (Ray et al., 2015). At higher
impact velocities, qualitatively different phenomena are observed. As shown in
figure 10, during the cavity expansion stage, the direction of the velocity
field is outward and upward around the enlarging cavity and the cavity
continues to expand even after the closure of the upper part. The cavity depth
reaches its maximum value at $t\approx 24$ ms as the cavity radius extends
continuously along its edge (see also figure 4e), which is very different from
the measurements for an even more energetic case impacting at $\sim 2U_{0}$ in
Engel (1966) where the crown and cavity arrive at their maximum positions at
about the same time. At the moment that the penetration depth reaches its
maximum value, the flow around the cavity bottom should have reached its
minimum value of essentially zero and shortly be redistributed by inertia
force, as seen at $t=24$ ms in figure 10. At approximately the same time
($t\approx 25$ ms), the downward spiral jet arrives at the cavity bottom and
penetrates deeply into the bulk, creating a subcavity that moves also
spirally, which may potentially accelerate the re-establishment of the
velocity field. By $t\approx 32$ ms, it can be observed that the velocity
field has been completely reversed in the pool and a new circulation has been
established. The flow around the cavity in the pool has transitioned from
outward-upward expansion to inward-downward contraction. Such a change of the
flow direction would result in a pressure buildup around the lower part of the
cavity, thus pushing the cavity to shallow. Air bubbles are entrapped in the
bulk in this process as can be seen at $t=32$, 38 and 48 ms in figure 10 and
figure 12(a). In the last frame of figure 10, an upward jet eventually rises
from the cavity floor and merges with the previous downward jet.

\begin{overpic}[scale={1.0}]{fig/Fig12.eps} \put(0.0,36.5){\small({a})}
\put(19.5,1.0){\scriptsize{$X$ $\rm{(mm)}$}}
\put(1.0,18.5){\scriptsize\rotatebox{90.0}{{$Y$ $\rm{(mm)}$}}}
\put(42.0,36.5){\small({b})} \put(51.7,36.3){\tiny{$U_{b}$}}
\put(51.7,34.5){\tiny{$Y_{b}$}} \put(68.0,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(42.0,21.0){\scriptsize\rotatebox{90.0}{{$Y_{b}/d$}}}
\put(97.0,21.0){\scriptsize\rotatebox{90.0}{{$U_{b}/U_{0}$}}} \end{overpic}
Figure 12: Kinematic behaviours of the entrapped large bubble. (a) Successive
positions of the vertical slices for the entrapped large bubble. The time
interval between each curve is 4 ms. (b) Time evolution of the vertical
centroid position (left) and the vertical speed (right) of the entrapped large
bubble. The bubble sinks at the expansion stage and then starts to shallow
from its bottom due to the concentric axial pressure, which eventually leads
to a floating air bubble above the pool surface.

Figure 12(a) draws the successive shapes of the entrapped large bubble and (b)
plots the temporal motion of the bubble centroid in the vertical direction,
showing the shallowing steps of the cavity. The bubble expands to a maximum
position and then contracts from its bottom, eventually generating a toroidal
air bubble floating at the top of the pool.

As for the final collapsing stage, the central column thickens and merges with
the outer bubble wall, creating a horseshoe-shaped bubble that transforms
later into a hemisphere. For the next long period of time (more than 300 ms),
this toroidal bubble will be stretched and thinned under the action of surface
tension and eventually ruptures due to instabilities. Subsequently, the film
cap recedes along the periphery of the “ruptured hole” from one side,
scattering fine scale liquid droplets from the “receding rim” to the air, as
shown in the experiments by Murphy et al. (2015) (figure 3). The production of
such tiny droplets from the receding films has been also studied by Lhuissier
& Villermaux (2012) and Dasouqi et al. (2021).

###  4.4 Energy budget

As reviewed by Veron (2015), one of the primary motivations for studying the
liquid droplets ejection is to estimate the exchange of momentum, heat and
eventually mass transfer through gas-liquid interface induced by the
interfacial processes such as large scale breaking wave event (Mostert et al.,
2022), small scale bubble bursting (Berny et al., 2022) and drop impact (Liang
& Mudawar, 2016) etc. In the present work, the kinetic energy $E_{k}$,
gravitational potential energy $E_{g}$ and surface potential energy $E_{s}$ of
the liquid pool are calculated as follows:

$E_{k}=\frac{1}{2}\displaystyle\int_{V_{p}}\rho\|\boldsymbol{U}\|^{2}dV_{p}$
(4) $E_{g}=\displaystyle\int_{V_{p}}\rho\boldsymbol{g}ydV_{p}-E_{g0}$ (5)
$E_{s}=\displaystyle\int_{S_{p}}\sigma dS_{p}-E_{s0}$ (6)

The integrals are computed over the volume ($V_{p}$) and surface ($S_{p}$) of
the largest liquid continuum in the domain, which means that the small
droplets detached from the pool and the gas phase are not included here. The
time point when the first contact occurs between drop and pool is defined as
$t=0$ ms. $E_{g0}$ and $E_{s0}$ are the gravitational and surface potential
energies of the liquid phase at $t=0$ ms. In addition, the kinetic,
gravitational and surface potential energies of the ejected droplets are
integrated over each droplet ($V_{p}^{i}$ and $S_{p}^{i}$, respectively) using
the same equations like 4, 5 and 6, and the sum of energy carried by small
droplets is represented as $E_{d}$. The total energy of the liquid phase is
therefore calculated as $E_{T}=E_{k}+E_{g}+E_{s}+E_{d}$. $E_{0}$ is the
initial kinetic energy introduced by the impacting drop.

\begin{overpic}[scale={1.3}]{fig/Fig13.pdf} \put(0.0,97.5){\small({a})}
\put(95.5,96.7){\tiny{$E_{k}$}} \put(95.5,94.7){\tiny{$E_{g}$}}
\put(95.5,92.4){\tiny{$E_{s}$}} \put(95.5,90.4){\tiny{$E_{d}$}}
\put(95.5,88.3){\tiny{$E_{T}$}}
\put(1.0,74.0){\scriptsize\rotatebox{90.0}{{$E/E_{0}$}}}
\put(0.0,50.0){\small({b})} \put(95.0,49.5){\tiny{$P_{y}$}}
\put(95.0,47.3){\tiny{$P_{xx}$}} \put(51.0,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(1.0,27.5){\scriptsize\rotatebox{90.0}{{$P/P_{0}$}}} \end{overpic} Figure
13: Energy budget for the process of high-speed drop impact onto deep liquid
pool. (a) Temporal evolution of the energy budget in liquid phase. The
energies are normalised by the initial kinetic energy of the impacting drop.
(b) Temporal evolution of the momentum components in the pool (largest liquid
continuum), normalised by the initial impact momentum. From left to right, the
vertical dotted lines indicate the time point for the closure of the upper rim
of the crown ($t\approx 12$ $\upmu$s), the connection of the downward spiral
jet to the cavity bottom ($t\approx 24$ $\upmu$s) and the formation of the
central upward jet from the cavity floor ($t\approx 40$ $\upmu$s).

Figure 13(a) plots the temporal evolution of the energy aspects throughout the
simulation. The calculated energies in the liquid phase are normalized by the
initial kinetic energy of the impacting drop $E_{0}$. Initially, at the time
the drop hits the pool, the system’s energy is mainly composed of the kinetic
energy of the impacting drop. As the drop moves downwards, a sharp decrease of
the kinetic energy is observed at the very early time of impact, which can be
associated with the immediate “prompt splash” of the liquid fragmentation
($E_{d}$) and the energy dissipation due to strong vortical entanglement near
the neck region. The gravitational potential energy remains insignificant at
this stage since the splashing appears mostly near the pool surface.
Subsequently, a period of rapid expansion of crown and cavity takes place (see
also figure 3 and 10), which is reflected by the almost linear decrease of the
kinetic energy and noticeable increase of the gravitational and surface
potential energies. The droplet energy $E_{d}$ reaches its maximum value at
approximately the same time when the droplet count reaches its peak during the
sustained droplet shedding stage (see also section 5). Afterwards, the number
of droplets reduces continuously and $E_{d}$ becomes insignificant. Here, the
maximum total droplet energy accounts for about $8\%$ of the initial kinetic
energy, which is quite close to the theoretical premises by Engel (1967) that
around $5\%$ of the initial impact energy is carried away by secondary
droplets. At the instant of closing (around $t\approx 12$ ms), more than
$60\%$ of the initial energy is still present in the liquid phase. Once the
crown necks in, changes of the slope are observed in both $E_{k}$ and $E_{s}$,
where the decreasing speed of $E_{k}$ slows down gradually and $E_{s}$ reaches
a plateau. $E_{g}$ reaches a maximum value at approximately the same time when
the cavity expands to the maximum position ($t\approx 24$ ms) and it decreases
hereafter due to the inverse flow around the bubble canopy (see figure 10).
Meanwhile, the kinetic energy in the pool keeps almost constant during the
initial retracting stage ($t\approx 24\sim 40$ ms) and increases slightly
later due to the protrusion of the broad upward jet from the cavity floor. The
total energy $E_{T}$ decreases monotonically throughout the entire process and
nearly half of the initial energy is eventually converted to other forms,
including viscous energy dissipation, airflow/bubble energy as well as the
energy contained by the removed droplets.

To better understand the kinematic behaviours of the liquid bulk, we also
extract the momentum components of the largest liquid continuum (except small
droplets) as shown in figure 13(b). The total vertical momentum $P_{y}$ can be
used to describe the overall vertical motion of the liquid bulk and the total
momentum on the right half of the bulk along the $x$-axis $P_{xx}$ can be used
to describe its horizontal expansion, which are both normalised by the initial
impact momentum of the drop $P_{0}$. For $t<25$ ms, a general outward-upward
motion of the liquid bulk is shown in the graph, where the maximum horizontal
and vertical expansion speeds are found at $t\approx 6$ ms and $t\approx 8$ ms
respectively. Starting from $t\approx 25$ ms, the bulk sinks downwards and the
liquid starts to flow concentrically towards the impact axis, squeezing the
large air bubble to float upwards as demonstrated in figure 12(b). A turning
point of $P_{y}$ is found at $t\approx 40$ ms, which indicates the onset of
the broad upward jet from the cavity floor.

###  4.5 Transportation of drop liquid

Figure 14: Transportation of the passive tracer for drop liquid. From left to
right and top to bottom, the corresponding times are 0.5, 3, 14, 20, 32 and 48
ms after impact. The scale bar is 10 mm long. Azimuthal destabilization is
captured at the edge of the drop liquid, which therefore produces secondary
droplets from its tips. At the later stage of impact, the thin drop film is
penetrated and mixed by the downward-moving jet, producing the drop liquid
plume deep inside the target pool.

In the simulation, a passive tracer field $f_{p}$ is initially added to the
impacting drop for the purpose of following the trace of the drop liquid.
Experimentally speaking, this can be achieved by adding colours/dye to the
liquids (Engel, 1966; Thoroddsen, 2002; Bisighini & Cossali, 2011). The tracer
field is then advected by the following equation:

$\frac{\p f_{p}}{\p t}+\boldsymbol{U_{f}}\bcdot\nabla f_{p}=0$ (7)

Figure 14 shows the distribution of drop liquid at different stages of impact.
Shortly after contact, the main part of the drop sits on the “pool cavity” and
the radially spreading thin film extends from its edges ($t=0.5$ ms). For the
next few milliseconds, the drop deforms and expands rapidly into a thin liquid
layer coated along the interior surface of the pool. Meanwhile, liquid threads
are emitted from the rim of the drop liquid and then break up into small
droplets on its tips ($t=3$ ms), which is comparable with the formation of
azimuthal destabilization on the retracting flattened drop edge at relatively
low impact energy reported in Lhuissier et al. (2013). By the time of closing,
it can be observed that these “drop threads” meet at the closure point and are
later transported backwards to the bulk by the downward-moving jet. At the
receding stage, the thin drop film starts to propagate towards the cavity
bottom due to the surface-tension capillary waves and is later transported
deep into the pool when the central spiral jet impinges the cavity bottom,
seen at $t=32$ ms and 48 ms in figure 14.

\begin{overpic}[width=433.62pt]{fig/Fig15.eps} \put(0.0,45.0){\small({a})}
\put(37.0,45.0){\small({b})} \put(71.0,44.3){\tiny{$T_{b}$(drop/pool
interface)}} \put(71.0,42.4){\tiny{$T_{a}$(air/drop interface)}}
\put(71.0,40.7){\tiny{$T_{\delta}$(drop thickness)}}
\put(71.0,38.8){\tiny{Analytical prediction \cite[citet]{\@@bibref{Authors
Phrase1YearPhrase2}{berberovic2009drop}{\@@citephrase{(}}{\@@citephrase{)}}})}}
\put(71.0,37.0){\tiny{Analytical prediction \cite[citet]{\@@bibref{Authors
Phrase1YearPhrase2}{bisighini2010crater}{\@@citephrase{(}}{\@@citephrase{)}}})}}
\put(69.0,2.0){\scriptsize{$tU_{0}/d$}}
\put(38.5,25.0){\scriptsize\rotatebox{90.0}{{$Y/d$}}} \end{overpic} Figure 15:
(a) Sketch of the drop penetration. Boundaries between different fluid
components are differentiated by isosurface $f_{p}=0.5$. The positions of the
upper point $T_{a}$, lower point $T_{b}$ and the thickness of the drop tracer
$T_{\delta}$ along the vertical axis of symmetry are tracked. (b) Time
variations of $T_{a}$, $T_{b}$ and $T_{\delta}$ along the axial direction. The
dashed line shows the asymptotic solution proposed by Berberović et al.
(2009). The solid line shows the theoretical estimation of the penetration
depth proposed by Bisighini et al. (2010). The dimensionless time $tU_{0}/d=2$
is indicated by vertical dotted line.

Different from the gas/liquid interface that can be recorded directly by the
camera, the “virtual” drop/pool interface is usually invisible if the fluids
in the drop and the receiving pool are the same, but the estimation of the
kinematic behaviours of this layer is vital for theoretical studies of crater
evolution (Bisighini & Cossali, 2011). As indicated in figure 15(a), the
boundaries between different components can be easily differentiated here by
the isosurface $f_{p}=0.5$. Figure 15(b) plots the temporal variation for the
positions of the upper point $T_{a}$, lower point $T_{b}$ as well as the
thickness of the drop tracer $T_{\delta}$ along the impact axis. As expected,
the upper point $T_{a}$ moves rapidly adjoining the air while the penetration
speed of the lower point $T_{b}$ is greatly decelerated by the reacting flows
from the target liquid, thus causing the decrease in drop thickness. It can be
seen that the drop deforms significantly and its thickness decreases
noticeably during an initial dimensionless time period $tU_{0}/d\leqslant 2$.
For the later stages, $T_{\delta}$ reaches a plateau and will not experience
much difference throughout the expansion stage. Slight fluctuations of
$T_{\delta}$ can be anticipated while the drop recedes into the cavity bottom
at the contraction stage.

Following previous investigations in literature, a critical dimensionless time
$tU_{0}/d\approx 2$ is usually used to subdivide the evolution of drop impact
into two phases (Fedorchenko & Wang, 2004). During the first phase
($tU_{0}/d\leqslant 2$), the drop deforms and extends above the bulk cavity,
and the interfaces of air/drop and drop/pool are clear-cut. The penetration
speed of the drop/target interface during this time can be approximated as
half of the impact speed $U_{p}\approx U_{0}/2$, which is well known from the
penetration mechanics (Birkhoff et al., 1948; Yarin et al., 1995) and has been
previously applied for the analytical study of drop impact (Fedorchenko &
Wang, 2004; Berberović et al., 2009). At times $tU_{0}/d>2$, the flow effects
in the thin drop layer are negligibly small and the cavity expansion thus can
be approximated as the shape of the drop/pool interface. Berberović et al.
(2009) developed a theoretical approach to estimate the penetration depth of
drop impact $T_{d}$ based only on the linear momentum balance of the liquid
around the cavity and gave an asymptotic solution as
$T_{d}=2^{-4/5}(5t-6)^{2/5}$ for $tU_{0}/d>2$ at high $Fr$, $We$ and $Re$
numbers. The predicted results using the above asymptotic formula are shown in
figure 15(b) (dashed line) and a fairly good agreement is found with the
simulated results. Since the effects of surface tension, viscosity and gravity
are neglected in this theory, it can not predict accurately the later stages
of impact where the deformation of the shape of the cavity becomes significant
due to gravity and capillary waves. Bisighini et al. (2010) proposed a
theoretical model based on the potential flow theory that accounted for the
effects of inertia, gravity, viscosity and surface tension for sufficiently
high Reynolds and Weber numbers, using the combination of the sphere expansion
and its translation along impact axis. A system of ordinary differential
equations is obtained and numerically solved using initial conditions:

$\ddot{\alpha}=-\frac{3}{2}\frac{\dot{\alpha}^{2}}{\alpha}-\frac{2}{\alpha^{2}We}-\frac{1}{Fr}\frac{\zeta}{\alpha}+\frac{7}{4}\frac{\dot{\zeta}^{2}}{\alpha}-\frac{4\dot{\alpha}}{\alpha^{2}Re}$
(8)
$\ddot{\zeta}=-3\frac{\dot{\alpha}\dot{\zeta}}{\alpha}-\frac{9}{2}\frac{\dot{\zeta}^{2}}{\alpha}-\frac{2}{Fr}-\frac{12\dot{\zeta}}{\alpha^{2}Re}$
(9)

Where $\alpha$ and $\zeta$ donate the dimensionless crater radius and axial
coordinate of the centre of the sphere and the dimensionless penetration depth
is expressed as $\alpha+\zeta$. As explained by Bisighini et al. (2010), the
initial conditions can be obtained from the initial phase ($tU_{0}/d\leqslant
2$) using the forms: $\dot{\alpha}\approx 0.17$,
$\alpha\approx\alpha+0.17\tau$, $\dot{\zeta}\approx 0.27$,
$\zeta\approx-\alpha_{0}+0.17\tau$, and the dimensionless width of the cavity
can be estimated using the geometrical conditions
$W=2\sqrt{\alpha^{2}-\zeta^{2}}\approx
2\sqrt{(\alpha_{0}+0.17\tau)^{2}-(0.27\tau-\alpha_{0})^{2}}$. By fitting the
simulated bulk cavity width using the least-mean-square method, the constant
is obtained as $\alpha_{0}=0.79$ in the present case, thus the initial
conditions for equations 8 and 9 are $\alpha(2)=1.13$, $\dot{\alpha}(2)=0.17$,
$\zeta(2)=-0.25$ and $\dot{\zeta}(2)=0.27$. As shown in the solid line in
figure 15(b), the temporal variation of the predicted depth of drop/pool
interface agrees very well with our numerical results during the expansion
stage of impact. As for the retraction stage, discrepancies between the
theoretical prediction and the simulation/experiment become more pronounced as
demonstrated in figure 4(e), which can be explained as the fact that the shape
of the cavity does not follow the spherical expansion anymore due to the
influence of the central spiral jet and the propagation of capillary waves,
the model therefore is invalid for the retraction phase.

##  5 Airborne Droplets

###  5.1 Source of secondary droplet

According to the investigation of Deegan et al. (2007) in the parameter range
$Re\leq 5000$ and $We\leq 1400$, at least three sources of secondary droplets
can be identified during the impact: (i) prompt instability of the ejecta
sheet occurring immediately after contact, which produces very small droplets,
(ii) rim instability of the ejecta sheet that produces medium-sized droplets
and (iii) rim instability of the crown that produces large droplets from
liquid jets. These different mechanisms are typically interdependent and the
earlier ones influence the later ones, which therefore further complicates the
characterization of these tiny airborne droplets.

\begin{overpic}[scale={0.75}]{fig/Fig16.pdf} \put(-3.0,98.5){\small({a})}
\put(-3.0,65.0){\small({b})} \put(-3.0,31.0){\small({c})} \end{overpic} Figure
16: Different mechanisms of droplet production at different stages of impact.
The images are shown under different magnifications. (a) $t=30$ $\upmu$s, the
“prompt splash” that occurs at the very early time of the impact near the neck
region due to irregular rupture and breakup of the ejecta film. (b) $t=650$
$\upmu$s, the sustained “crown splash” due to breakup of the thin liquid
ligaments on the top of the crown rim. (c) Partially resolved tiny droplets
near the pool surface produced by secondary impact/bubble bursting. The left
panel shows the locations of splashing and the right panel demonstrates the
mesh structures nearby.

Figure 16 qualitatively illustrates some primary mechanisms of droplet
production during the impact captured by simulation. The left panel shows the
stages and locations where tiny droplets are generated and the right panel
demonstrates the mesh structures on the overlaid section. As shown in figure
16(a), a great number of microdroplets are generated near the neck region
immediately after the first contact due to the very early breakups of ejecta
film, so-called the “prompt splash” (Deegan et al., 2007; Marcotte et al.,
2019). We have estimated that around $85\%$ of droplets produced at this time
step have the equivalent mean diameter of more than double the size of the
smallest cell and $70\%$ of them are in the size range of $10\sim 30$
$\upmu$m, suggesting that the statistical analysis of the initial smallest
droplets should be treated carefully. After the intricate early splashing, a
“smooth” liquid sheet rises around the contact line to form the cylindrical
crown, emanating thin liquid ligaments on its rims from a fairly regular
distance. A sustained droplet ejection, so-called the “crown splash”, is
therefore observed. The size of droplets produced at this stage is greatly
influenced by the thickness of ligaments and increases with time. Figure 16(b)
shows the fragmentation on the tips of the ligaments at $t=650$ $\upmu$s. From
the right panel, it can be clearly seen that the droplets produced at the
present time instant are sufficiently larger than the minimum scale of the
cell. Figure 16(c) shows the production of smaller droplets from the
“secondary impact” caused by the previous generation. When the first-born
droplets fall back and impinge to the pool, it may produce another splashing
or breakups, which are generally partially resolved. The radii of these
smallest droplets are represented approximately by the smallest size of the
mesh (right panel). Besides the above sources, very small child-drops can be
also ejected from bubble bursting, which are usually not fully resolved as
same as figure 16(c). Fewer larger droplets are also observed later from the
downward-moving central jet after the closure of the upper crown.

###  5.2 Droplet statistics

\begin{overpic}[width=433.62pt]{fig/Fig17.pdf} \put(0.0,34.7){\small({a})}
\put(25.5,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(1.0,20.0){\scriptsize\rotatebox{90.0}{{$N_{d}$}}}
\put(49.5,34.7){\small({b})} \put(75.5,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(50.5,18.0){\scriptsize\rotatebox{90.0}{{$M_{d}/M_{0}$}}} \end{overpic}
Figure 17: (a) Temporal evolution of the total number of secondary droplets.
The vertical dotted line indicates the maximum droplet count at $t\approx 650$
$\upmu$s. (b) Temporal evolution of the total mass of secondary droplets. The
time point when the impacting drop touches the pool is chosen as the reference
time. The total mass of secondary droplets $M_{d}$ is scaled by the mass of
the impacting drop $M_{0}$.

We now discuss the statistics of droplets. Here we only analyze the droplet
information captured with higher resolutions ($L_{max}$=13 and 14) for the
first 4 ms after impact. The time variation of the total number of droplets is
plotted in figure 17(a). It should be noted that a short time disturbance is
presented during $t\approx 0.21\sim 0.27$ ms on the curve, which might be
explained as the loss of the smallest droplets at $L_{max}$=14 within the
extra refinement layer near the surface of the pool (see section 2.3), as it
happens at approximately the same time when many droplets reach this height
and the extra refinement layer is removed. Neglecting the artefacts at the
early time, a unimodal curve of droplet production is found. Shortly after the
first contact, a large number of very fine droplets are scattered from the
ruptured liquid ejecta (corresponding to figure 16a), which is reflected in
figure 17(a) for the sharp increase of the droplet count at the beginning.
Subsequently, larger droplets are emitted almost continuously from the tips of
the ligaments. The maximum droplet count is reached 650 $\upmu$s after impact
with around 4340 droplet population, and the dynamics around this time are
qualitatively shown in figure 16(b). As the crown grows, ligaments merge along
the rim and become shorter and thicker, producing larger droplets till the
closure of the upper part. A change of the decreasing slope can be found at
$t\approx$ 1.7 ms in figure 17(a), which indicates the timing when lots of
droplets start to exit the field of view and are removed from the
computational domain.

Figure 17(b) plots the time evolution of the mass ratio of the total secondary
droplets ($M_{d}$) to the impact drop ($M_{0}$). An overall increasing trend
of the mass transfer from pool to air is found. The number of droplets starts
to decrease after the peak ($t\approx 650$ $\upmu$s) but the total mass of
droplets keeps increasing, which also reveals the fact that the droplets
produced at the later “crown splash” stage are in much larger scales than the
early splashing.

\begin{overpic}[width=433.62pt]{fig/Fig18.eps} \put(0.0,41.0){\small({a})}
\put(25.0,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(1.5,22.0){\scriptsize\rotatebox{90.0}{{$S_{d}$ $\rm{(\upmu m)}$}}}
\put(54.5,41.0){\small({b})} \put(87.0,40.3){\tiny{t = 30 $\rm{\upmu s}$}}
\put(87.0,38.2){\tiny{t = 150 $\rm{\upmu s}$}} \put(87.0,36.0){\tiny{t = 650
$\rm{\upmu s}$}} \put(87.0,33.8){\tiny{$\overline{t}=3\sim 4$ $\rm{ms}$}}
\put(78.0,1.0){\scriptsize{{$\overline{N}(S_{d},t)$}}}
\put(55.5,22.0){\scriptsize\rotatebox{90.0}{{$S_{d}$ $\rm{(\upmu m)}$}}}
\end{overpic} Figure 18: (a) Temporal contour of droplet size distribution
during $0\sim 4$ ms of impact. The horizontal dashed lines (red) indicate the
length scale of doubling the minimum cell size $\sim 2\Delta$. (b) Droplet
size distribution at different time slices in (a). Here it shows $t=30$
$\upmu$s and $t=150$ $\upmu$s when the liquid ejecta ruptures and produces
very tiny droplets; $t=650$ $\upmu$s when the droplet count reaches the
maximum. Large size droplets ($S_{d}>100$ $\upmu$m) are only generated from
the ligament breakups at the “crown splash” stage. Time-averaged droplet size
distribution is calculated using 9 time slices from the time window $t=3\sim
4$ ms. A bi-model distribution of droplet size is found.

Figure 18(a) plots the contour map of droplet size distribution per bin size
$\Delta r$ at various times. The equivalent mean diameter of each droplet
$S_{d}$ is calculated as a sphere using the integrated volume fraction of the
liquid phase. Initially, a lot of relatively small droplets, $S_{d}$ ranging
around $8\sim 25$ $\upmu$m, are generated immediately upon impact due to the
mechanism of “prompt splash” (see figure 16a). In the experiments, a ring of
fine spray containing mostly $6\sim 19$ $\upmu$m droplets was also recorded at
the instant of impact, which agrees well with our numerical results, but the
existence of smaller droplets was not able to be determined owing to limited
camera resolution. In the simulation, a sharp increase of much smaller
droplets ($S_{d}\leqslant 8$ $\upmu$m) is captured at the initial splashing
stage ($t<200$ $\upmu$s), which can be associated with the partially resolved
very small droplets from the rupture and breakup of the early ejecta film and
“secondary impact”. It is worth noting that the minimum cell scales are
$\Delta\approx 3.9$ $\upmu$m at stage $S1$ and $\Delta\approx 7.8$ $\upmu$m at
stage $S2$ in the present simulation as introduced in 2.3, and we assume that
the geometries can not be represented for those droplets whose radii are
smaller than the minimum cell size (red dashed line in figure 18a), which
could influence the breakup physics at smaller scales. However, grid
convergence of the numerical atomization remains an open question as has been
discussed by Herrmann (2013) and Pairetti et al. (2020), although the latter
author suggested that 8 cells per drop diameter would probably be enough to
preserve most of their physical behaviours. Referring to the recent works of
Wang et al. (2016) and Mostert et al. (2022) concerning droplet spray in the
action of breaking waves, droplets radii greater than $2\Delta$ would be able
to maintain a spherical shape and are approximately grid-converged with proper
time-averaging procedures, suggesting that at least $4\sim 8$ points per
droplet diameter are essential to obtain a numerical convergence.

After the initial splashing, an “intact” liquid sheet rises cylindrically from
the contact line and disseminates secondary droplets from the thin liquid
ligaments as shown in figure 16(b). The size of droplets increases in time
with the gradually thickened ligaments at this period. In figure 18(a), three
main size concentrations become pronounced 2 ms after impact, which can be
briefly classified as the small-sized class $S_{d}\leqslant 16$ $\upmu$m, the
medium-sized class $16<S_{d}<50$ $\upmu$m and the large-sized class
$S_{d}\geqslant 100$ $\upmu$m. The first concentration range is associated
with the partially resolved smallest droplets and their count decreases
noticeably in time, indicating that they are mainly generated at the very
beginning of the impact and lots of them fall back to the bulk very quickly
after birth. The second concentration range lies in $16\sim 50$ $\upmu$m
(moderate size) and they are the ones produced the most throughout the entire
simulation. The last concentration range is composed of large droplets (larger
than 100 $\upmu$m), increasing and narrowing in time, which follows the
experimental observation (Murphy et al., 2015) that fewer but larger droplets
will be detached from the gradually thickened ligaments with time.

Various Instantaneous time slices of droplet size distribution are provided in
figure 18(b). At $t\approx 150$ $\upmu$s, the “bumping” event between the
ejecta and the drop/pool just stopped and the primary mechanism of droplet
production starts to transition from “prompt splash” to “crown splash”.
Comparing the size distribution at $t=150$ $\upmu$s and $t=650$ $\upmu$s, very
similar profiles can be found in the group of moderate droplets ($16\sim 50$
$\upmu$m), which confirms that the medium-sized concentration in figure 18(a)
is composed primarily of the droplets from the mechanism (ii) summarized in
Deegan et al. (2007): early rim instability from the very thin ejecta sheet.
The discrepancy in small-sized droplets between these two time steps is
probably caused by the change of mesh resolution, since the smallest cell size
at $t>200$ $\upmu$s has been doubled. Large-sized droplets of $S_{d}>100$
$\upmu$m appear at $t>150$ $\upmu$s, suggesting that large scale droplets are
only produced by the fragmentation of thin ligaments at the “crown splash”
stage.

As the ligaments merge and thicken, the newly detached droplets increase
gradually in size, and the most abundant droplet information occurs around
$3\sim 4$ ms after impact (Murphy et al., 2015). The time-averaged droplet
size distribution over the time window $t=3\sim 4$ ms is therefore calculated.
As shown in figure 18(b), a bimodal distribution of droplet size can be
observed from the histogram, corroborating the separated size concentration
ranges in figure 18(a). The bimodal size distribution of liquid sheet
fragmentation has been also reported by some other studies (Afeti & Resch,
1990; Villermaux & Bossa, 2009, 2011; Murphy et al., 2015). Two primary
distribution ranges can be briefly distinguished. The first one is composed of
small-sized and medium-sized droplets peaking around 30 $\upmu$m. The second
primary range appears for the large scale droplets ranging from 100 $\upmu$m
to 500 $\upmu$m.

\begin{overpic}[width=433.62pt]{fig/Fig19.eps} \put(0.0,41.0){\small({a})}
\put(25.0,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(1.5,22.0){\scriptsize\rotatebox{90.0}{{$S_{d}$ $\rm{(\upmu m)}$}}}
\put(54.0,41.0){\small({b})} \put(78.0,1.0){\scriptsize{$t$ $\rm{(ms)}$}}
\put(55.5,23.0){\scriptsize\rotatebox{90.0}{{$N_{d}$}}} \end{overpic} Figure
19: Statistics of droplets who tend to re-merge with the liquid bulk. (a)
Temporal contour of the droplet size distribution for the “re-merging”
droplets. (b) Evolution of the “re-merging” droplets count with time.

Lastly, we would like to focus on those droplets that tend to fall back to the
liquid bulk, which is vital for the estimation of mass transfer through the
air-sea interface. In practice, those droplets that move along the impact
direction (downwards) and whose centroid are located less than 100 $\upmu$m to
the free surface of the liquid pool are selected and assumed as the ones that
will re-join the pool. Figure 19(a) shows the temporal contour of the size
distribution for the “re-merging” droplets and (b) shows its count in time. It
can be found that most of the initially produced small-sized droplets from
“prompt splash” vibrate near the free surface and tend to re-join the bulk
shortly after generation. At the beginning of the “crown splash”, very thin
ligaments stretch out from the rims of the downward bending ejecta sheet and
send out tiny droplets towards the target liquid as demonstrated in figure
3(b), which thereby contribute to a second primary peak of the medium-sized
droplets in figure 19(a,b). As the impact advances, the direction of the
liquid ligaments transients quickly from horizontal-outward to vertical-upward
and the droplets are generally pinched off upwards in larger angles
accordingly. The amount of the “re-merging” droplet becomes insignificant 1 ms
after impact. Finally, it is worth noting that the presence of wind may
significantly advect the smallest droplets and delay or prevent them from
rejoining the bulk.

##  6 Conclusion

In this work, the high-energy splash of drop impact onto a deep volume of the
same liquid pool has been investigated with high-resolution direct numerical
simulation (DNS) in 3D. The calculations have been conducted in the exact same
configurations studied previously by Murphy et al. (2015), in order to perform
detailed comparisons with the experiments and prepare an in-depth analysis of
the splashing dynamics. The BASILISK open-source solver, which combines a VOF
description of the gas/liquid interfaces and an adaptive octree grid
refinement, has been employed to simulate the process of drop impact.
Qualitative and quantitative comparisons between numerical and experimental
results have been conducted in terms of the morphological behaviours of
splashing, kinematics of crown and cavity as well as the distributions of
secondary droplets through a particular field of view, which efficiently
validate the present numerical strategies and therefore enable the discussion
of the internal mechanisms for high-speed drop impact afterwards.

Following the experimental observations of Murphy et al. (2015), we performed
a detailed investigation on flow physics and splashing behaviours of high-
speed drop impact, serving as an important supplementary study for this issue.
Firstly, the very early instantaneous motions in the vicinity of the neck
region were discussed under sufficient time resolution. We have confirmed the
existence of two different mechanisms of air entrapment in the neck of
connection, namely the entrapment of axisymmetric bubble rings driven by high
localized pressure at $R_{n}/R_{d}<35\%$ (Oguz & Prosperetti, 1989) and the
entrapment of isolated bubbles/bubble rings due to the unstable oscillations
of the ejecta base (Weiss & Yarin, 1999). Moreover, the calculation
successfully captured the initialization of the irregular azimuthal
undulations along the outer edge of the neck, which thereby breaks the
axisymmetry of the motions in microscopic scales. Thereafter, detailed
information on the internal flows such as velocity, pressure as well as energy
budget was extracted from the calculation, to further explain the
corresponding physical phenomena observed in experiments. We showed that
azimuthal destabilization occurs on the edges of the flattening drop, breaking
up on its tips and participating in the production of child-droplets. These
“liquid threads”, emitted from the initial impact drop, grow together with the
uprising crown and meet upon closure, and they are eventually transported
backwards to the bulk along with the penetration of the central spiral jet.
Lastly, we presented the statistics of airborne droplets produced by the high-
energy drop-pool collision. The results showed that a great number of
microdroplets are produced immediately after contact by the irregular “prompt
splash” within the time scale of $t<200$ $\upmu$s, composing the most
populated small and moderate sizes. The earliest tiny droplets are sprayed
right above the surface of the pool and most of them move towards the target
liquid, suggesting the great possibility that a big part of them may return to
the liquid bulk shortly after birth. Large droplets ($S_{d}>100$ $\upmu$m) are
only observed from the fragmentation of the rim ligaments during the “crown
splash” stage, which therefore forms a gradually increasing but narrowing
distribution ribbon, reflecting the merging and thickening activities of the
ligaments on the rim of the crown. Finally, the bimodal size distribution of
secondary droplets has been found.

We are aware that the splashing dynamics induced by the most energetic regime
of drop impact are way more complicated than what has been mentioned here.
Further analysis of the multi-scale flow physics needs to be conducted under
sufficient spatial and temporal resolutions in the future. Details of the
bubble ring entrapment are still not well understood, as well as its motions
and breakups. The physical mechanism that is responsible for the early
azimuthal instability remains unknown. Recent experimental observations
(Thoraval et al., 2013; Li et al., 2018) have suggested that these early
dynamics may be greatly affected by the intricate vortical motions and three-
dimensional instabilities, which certainly added more complexities to the
analysis. For drop impact on the same liquid pool, large parameter space of
$Re$ and $We$ need to be studied to determine the critical conditions of the
large bubble entrapment (bubble canopy regime), as well as to explore its
influence on the closure event (closure time, closure height, large bubble
volume, floating bubble and its burst). The formation and motion of the
central spiral jet also need to be analyzed in detail and compared with
various types of liquid column jets. Statistics of droplets and bubbles need
to be collected under various impacting conditions to form a more
comprehensive database, which could therefore reflect directly the
mass/momentum transfer through gas-liquid interfaces as well as enlighten
various applications where this process is involved (oil dispersant, spray
cooling, metallurgy, etc).

[Acknowledgements]We appreciate beneficial discussions and help from the
BASILISK community. Simulations were performed using computational resources
on Advanced Research Computing (ARC) at Virginia Tech.

[Funding]This work is supported by the scholarship from China Scholarship
Council (CSC) under the Grant CSC NO. 201908320462.

[Declaration of interests]The authors report no conflict of interest.

[Author ORCID]Hui Wang, https://orcid.org/0000-0001-9733-0150; Shuo Liu,
https://orcid.org/0000-0002-8530-8359.

##  Appendix A Effects of spatial resolution

\begin{overpic}[scale={1.5}]{fig/Fig20.pdf} \put(0.0,98.4){\small({a})}
\put(0.0,86.2){\small({b})} \put(0.0,74.3){\small({c})}
\put(13.0,44.0){\small({d})} \put(64.0,43.7){\tiny{$L_{max}$ = $11$}}
\put(64.0,41.5){\tiny{$L_{max}$ = $12$}} \put(64.0,39.4){\tiny{$L_{max}$ =
$13$}} \put(64.0,37.2){\tiny{$L_{max}$ = $14$}} \put(43.0,2.0){\scriptsize{$t$
$(\upmu\rm{s})$}} \put(14.0,25.0){\scriptsize\rotatebox{90.0}{{$N_{c}$}}}
\end{overpic} Figure 20: Effects of the maximum mesh refinement level on
splashing behaviours. (a) Evolution of early-time splashing obtained with
$L_{max}=12$, showing 10, 40, 70 and 100 $\upmu$s after impact. (b) Evolution
of early-time splashing obtained with $L_{max}=14$, showing 10, 40, 70 and 100
$\upmu$s after impact. (c) Comparison of air-water interfaces at $t=50$
$\upmu$s captured by calculation with $L_{max}=12$ (left), experiment (middle)
and calculation with $L_{max}=14$ (right). (d) Time evolution of cell numbers
under different maximum mesh refinement levels.

In this appendix, we discuss the effects of the minimum spatial resolution.
Using the Adaptive Mesh Refinement (AMR) algorithm, preliminary tests have
been carried out by varying the maximum refinement level $L_{max}$ from 11 to
14, corresponding to the minimum cell size of $3.9\sim 31.25$ $\upmu$m in
reality. Two types of distinctive splashing phenomena are observed under
different spatial resolutions. For simulations at lower levels ($L_{max}$=12
in figure 20a), a thin layer of ejecta sheet rises smoothly from the contact
line and later becomes the leading edge of the crown, and secondary droplets
detach nearly axisymmetrically from the rim of the liquid sheet. Increasing
the maximum mesh refinement to higher levels ($L_{max}$=14 in figure 20b), a
disturbed incoherent liquid sheet accompanied by more “randomly” distributed
secondary droplets are captured (details in section 4.1 and 5). Referring the
splashing classifications in Thoroddsen et al. (2011) and Thoraval et al.
(2012) based on the dimensionless Ohnesorge number ($Oh=\mu/\rho\sigma d$) and
Splash number ($K=We\sqrt{Re}$), irregular splashing always occurs for the
most energetic cases. As the parameters in our present case are much higher
than their range of study, thinner liquid sheet and more complex interfacial
deformations/breakups can be even expected.

Figure 20(c) compares the shape of the air-water interface captured with
$L_{max}$=12 (left), high-speed camera (middle) and $L_{max}$=14 (right) 50
$\upmu$s after impact. The emerge-ruptured liquid sheet together with a great
number of irregularly distributed tiny droplets calculated at $L_{max}$=14 are
consistent well with the experimental observations, while a rather smooth
ejecta is presented with $L_{max}$=12 when the spatial resolution is
insufficient. As a formal study of mesh influence could not be performed,
given the computational resources required, the successful reproduction of the
primary features of the early splashing confirms that a maximum refinement
level at $L_{max}=14$ is essential for capturing the “correct” physical
dynamics (see also section 3.1).

It should be noted that each time the maximum refinement level is increased,
eight “children” cells will be divided from the tree-based structure in 3D.
Although the adaptive wavelet algorithm concentrates the smallest scales
mostly near the interfaces of the impacting area, the computational
requirements remain harsh at high levels. Figure 20(d) shows the time
evolution of the cell number for the first 250 $\upmu$s of the simulations
under different spatial resolutions. It can be observed that the total number
of cells is more than double with each additional increase of the $L_{max}$,
and this number will become considerably huge later when the large crater and
the substantial droplets and bubbles are developed.

## References

  * Afeti & Resch (1990) Afeti, GM & Resch, FJ 1990 Distribution of the liquid aerosol produced from bursting bubbles in sea and distilled water. Tellus B 42 (4), 378–384.
  * Agbaglah et al. (2015) Agbaglah, Gilou, Thoraval, Marie-Jean, Thoroddsen, Sigurdur T, Zhang, Li V, Fezzaa, Kamel & Deegan, Robert D 2015 Drop impact into a deep pool: vortex shedding and jet formation. J. Fluid Mech 764, R1.
  * Bell et al. (1989) Bell, John B, Colella, Phillip & Glaz, Harland M 1989 A second-order projection method for the incompressible navier-stokes equations. Journal of Computational Physics 85 (2), 257–283.
  * Berberović et al. (2009) Berberović, Edin, van Hinsberg, Nils P, Jakirlić, Suad, Roisman, Ilia V & Tropea, Cameron 2009 Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Physical review E 79 (3), 036306.
  * Berny et al. (2022) Berny, Alexis, Deike, Luc, Popinet, Stéphane & Séon, Thomas 2022 Size and speed of jet drops are robust to initial perturbations. Physical Review Fluids 7 (1), 013602\.
  * Birkhoff et al. (1948) Birkhoff, Garrett, MacDougall, Duncan P, Pugh, Emerson M & Taylor, Sir Geoffrey 1948 Explosives with lined cavities. Journal of Applied Physics 19 (6), 563–582.
  * Bisighini & Cossali (2011) Bisighini, Alfio & Cossali, Gianpietro Elvio 2011 High-speed visualization of interface phenomena: single and double drop impacts onto a deep liquid layer. Journal of visualization 14 (2), 103–110.
  * Bisighini et al. (2010) Bisighini, Alfio, Cossali, Gianpietro E, Tropea, Cameron & Roisman, Ilia V 2010 Crater evolution after the impact of a drop onto a semi-infinite liquid target. Physical Review E 82 (3), 036319.
  * Boulton-Stone & Blake (1993) Boulton-Stone, JM & Blake, JR 1993 Gas bubbles bursting at a free surface. Journal of Fluid Mechanics 254, 437–466.
  * Brackbill et al. (1992) Brackbill, Jeremiah U, Kothe, Douglas B & Zemach, Charles 1992 A continuum method for modeling surface tension. Journal of computational physics 100 (2), 335–354.
  * Brambilla & Guardone (2013) Brambilla, Paola & Guardone, Alberto 2013 Automatic tracking of corona propagation in three-dimensional simulations of non-normal drop impact on a liquid film. Computing 95 (5), 415–424.
  * Brambilla & Guardone (2015) Brambilla, Paola & Guardone, Alberto 2015 Assessment of dynamic adaptive grids in volume-of-fluid simulations of oblique drop impacts onto liquid films. Journal of Computational and Applied Mathematics 281, 277–283.
  * Castillo Orozco (2015) Castillo Orozco, Eduardo 2015 Droplet impact on liquid pools: secondary droplets formation from rayleigh jet break-up and crown splash .
  * Castillo-Orozco et al. (2016) Castillo-Orozco, Eduardo, Davanlou, Ashkan, Choudhury, Pretam K & Kumar, Ranganathan 2016 On the impact of liquid drops on immiscible liquids. In ASME 2016 14th International Conference on Nanochannels, Microchannels, and Minichannels collocated with the ASME 2016 Heat Transfer Summer Conference and the ASME 2016 Fluids Engineering Division Summer Meeting. American Society of Mechanical Engineers Digital Collection.
  * Castrejón-Pita et al. (2012) Castrejón-Pita, Alfonso A, Castrejón-Pita, JR & Hutchings, Ian M 2012 Experimental observation of von kármán vortices during drop impact. Physical Review E 86 (4), 045301.
  * Chandrasekhar (2013) Chandrasekhar, Subrahmanyan 2013 Hydrodynamic and hydromagnetic stability. Courier Corporation.
  * Cheng & Lou (2015) Cheng, M & Lou, J 2015 A numerical study on splash of oblique drop impact on wet walls. Computers & Fluids 115, 11–24.
  * Dasouqi et al. (2021) Dasouqi, Ali A, Yeom, Geum-Su & Murphy, David W 2021 Bursting bubbles and the formation of gas jets and vortex rings. Experiments in Fluids 62 (1), 1–18.
  * Deegan et al. (2007) Deegan, RD, Brunet, P & Eggers, J 2007 Complexities of splashing. Nonlinearity 21 (1), C1.
  * Deegan et al. (2008) Deegan, Robert D, Brunet, Philippe & Eggers, Jens 2008 Rayleigh-plateau instability causes the crown splash. arXiv preprint arXiv:0806.3050 .
  * Deka et al. (2017) Deka, Hiranya, Ray, Bahni, Biswas, Gautam, Dalal, Amaresh, Tsai, Pei-Hsun & Wang, An-Bang 2017 The regime of large bubble entrapment during a single drop impact on a liquid pool. Physics of Fluids 29 (9), 092101.
  * Deng et al. (2007) Deng, Q, Anilkumar, AV & Wang, TG 2007 The role of viscosity and surface tension in bubble entrapment during drop impact onto a deep liquid pool. Journal of Fluid Mechanics 578, 119.
  * Duchemin et al. (2003) Duchemin, L, Eggers, Jens & Josserand, C 2003 Inviscid coalescence of drops. Journal of Fluid Mechanics 487, 167–178.
  * Eggers et al. (1999) Eggers, Jens, Lister, John R & Stone, Howard A 1999 Coalescence of liquid drops. Journal of Fluid Mechanics 401, 293–310.
  * Engel (1966) Engel, Olive G 1966 Crater depth in fluid impacts. Journal of Applied Physics 37 (4), 1798–1808.
  * Engel (1967) Engel, Olive G 1967 Initial pressure, initial flow velocity, and the time dependence of crater depth in fluid impacts. Journal of applied physics 38 (10), 3935–3940.
  * Ervik et al. (2014) Ervik, Åsmund, Hellesø, Svein Magne, Munkejord, Svend Tollak & Müller, Bernhard 2014 Experimental and computational studies of water drops falling through model oil with surfactant and subjected to an electric field. In 2014 IEEE 18th International Conference on Dielectric Liquids (ICDL), pp. 1–6. IEEE.
  * Fedorchenko & Wang (2004) Fedorchenko, Alexander I & Wang, An-Bang 2004 On some common features of drop impact on liquid surfaces. Physics of Fluids 16 (5), 1349–1365.
  * Ferreira et al. (1985) Ferreira, Alfredo G, Larock, Bruce E & Singer, Michael J 1985 Computer simulation of water drop impact in a 9.6-mm deep pool. Soil Science Society of America Journal 49 (6), 1502–1507.
  * Francois et al. (2006) Francois, Marianne M, Cummins, Sharen J, Dendy, Edward D, Kothe, Douglas B, Sicilian, James M & Williams, Matthew W 2006 A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. Journal of Computational Physics 213 (1), 141–173.
  * Fudge et al. (2021) Fudge, Ben D, Cimpeanu, Radu & Castrejón-Pita, Alfonso A 2021 Dipping into a new pool: The interface dynamics of drops impacting onto a different liquid. Physical Review E 104 (6), 065102.
  * Fuster & Popinet (2018) Fuster, Daniel & Popinet, Stéphane 2018 An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension. Journal of Computational Physics 374, 752–768.
  * Gekle & Gordillo (2010) Gekle, Stephan & Gordillo, José Manuel 2010 Generation and breakup of worthington jets after cavity collapse. part 1. jet formation. Journal of fluid mechanics 663, 293–330.
  * Gielen et al. (2017) Gielen, Marise V, Sleutel, Pascal, Benschop, Jos, Riepen, Michel, Voronina, Victoria, Visser, Claas Willem, Lohse, Detlef, Snoeijer, Jacco H, Versluis, Michel & Gelderblom, Hanneke 2017 Oblique drop impact onto a deep liquid pool. Physical review fluids 2 (8), 083602.
  * Guildenbecher et al. (2014) Guildenbecher, Daniel R, Engvall, Luke, Gao, Jian, Grasser, Thomas W, Reu, Phillip L & Chen, Jun 2014 Digital in-line holography to quantify secondary droplets from the impact of a single drop on a thin film. Experiments in fluids 55 (3), 1–9.
  * Harlow & Shannon (1967) Harlow, Francis H & Shannon, John P 1967 The splash of a liquid drop. Journal of Applied Physics 38 (10), 3855–3866.
  * Hendrix et al. (2016) Hendrix, Maurice HW, Bouwhuis, Wilco, van der Meer, Devaraj, Lohse, Detlef & Snoeijer, Jacco H 2016 Universal mechanism for air entrainment during liquid impact. Journal of fluid mechanics 789, 708–725.
  * Herrmann (2013) Herrmann, Marcus 2013 On simulating primary atomization. Atomization and Sprays 23 (12).
  * Hogrefe et al. (1998) Hogrefe, J Errett, Peffley, Nicholas L, Goodridge, Christopher L, Shi, William T, Hentschel, H George E & Lathrop, Daniel P 1998 Power-law singularities in gravity-capillary waves. Physica D: Nonlinear Phenomena 123 (1-4), 183–205.
  * van Hooft et al. (2018) van Hooft, J Antoon, Popinet, Stéphane, van Heerwaarden, Chiel C, van der Linden, Steven JA, de Roode, Stephan R & van de Wiel, Bas JH 2018 Towards adaptive grids for atmospheric boundary-layer simulations. Boundary-layer meteorology 167 (3), 421–443.
  * Howison et al. (2005) Howison, SD, Ockendon, JR, Oliver, JM, Purvis, R & Smith, FT 2005 Droplet impact on a thin fluid layer. Journal of Fluid Mechanics 542, 1–23.
  * Jain et al. (2019) Jain, Utkarsh, Jalaal, Maziyar, Lohse, Detlef & van der Meer, Devaraj 2019 Deep pool water-impacts of viscous oil droplets. Soft matter 15 (23), 4629–4638.
  * Jamali et al. (2020) Jamali, M, Rostamijavanani, A, Nouri, NM & Navidbakhsh, M 2020 An experimental study of cavity and worthington jet formations caused by a falling sphere into an oil film on water. Applied Ocean Research 102, 102319.
  * Jian et al. (2020) Jian, Zhen, Channa, Murad Ali, Kherbeche, Abderrahmane, Chizari, Hossain, Thoroddsen, Sigurdur T & Thoraval, Marie-Jean 2020 To split or not to split: dynamics of an air disk formed under a drop impacting on a pool. Physical Review Letters 124 (18), 184501.
  * Josserand & Zaleski (2003) Josserand, Christophe & Zaleski, Stephane 2003 Droplet splashing on a thin liquid film. Physics of fluids 15 (6), 1650–1657.
  * Khokhlov (1998) Khokhlov, Alexei M 1998 Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations. Journal of Computational Physics 143 (2), 519–543.
  * Kim et al. (2021) Kim, Dohyung, Lee, Jinseok, Bose, Arijit, Kim, Ildoo & Lee, Jinkee 2021 The impact of an oil droplet on an oil layer on water. Journal of Fluid Mechanics 906.
  * Leng (2001) Leng, Liow Jong 2001 Splash formation by spherical drops. Journal of Fluid Mechanics 427, 73–105.
  * Lherm et al. (2021) Lherm, Victor, Deguen, Renaud, Alboussière, Thierry & Landeau, Maylis 2021 Rayleigh-taylor instability in drop impact experiments. Physical Review Fluids 6 (11), 110501\.
  * Lhuissier et al. (2013) Lhuissier, Henri, Sun, Chao, Prosperetti, Andrea & Lohse, Detlef 2013 Drop fragmentation at impact onto a bath of an immiscible liquid. Physical review letters 110 (26), 264503\.
  * Lhuissier & Villermaux (2012) Lhuissier, Henri & Villermaux, Emmanuel 2012 Bursting bubble aerosols. Journal of Fluid Mechanics 696, 5–44.
  * Li et al. (2018) Li, EQ, Thoraval, Marie-Jean, Marston, JO & Thoroddsen, Sigurdur T 2018 Early azimuthal instability during drop impact. Journal of Fluid Mechanics 848, 821–835.
  * Li et al. (2019) Li, Jianxin, Zhang, Huang & Liu, Qianfeng 2019 Characteristics of secondary droplets produced by a single drop impacting on a static liquid film. International Journal of Multiphase Flow 119, 42–55.
  * Liang & Mudawar (2016) Liang, Gangtao & Mudawar, Issam 2016 Review of mass and momentum interactions during drop impact on a liquid film. International Journal of Heat and Mass Transfer 101, 577–599.
  * Liu (2018) Liu, Xinan 2018 Experimental study of drop impact on deep-water surface in the presence of wind. Journal of Physical Oceanography 48 (2), 329–341.
  * Macklin & Hobbs (1969) Macklin, WC & Hobbs, PV 1969 Subsurface phenomena and the splashing of drops on shallow liquids. Science 166 (3901), 107–108.
  * Manzello & Yang (2002) Manzello, Samuel L & Yang, Jiann C 2002 An experimental study of a water droplet impinging on a liquid surface. Experiments in fluids 32 (5), 580–589.
  * Marcotte et al. (2019) Marcotte, Florence, Michon, Guy-Jean, Séon, Thomas & Josserand, Christophe 2019 Ejecta, corolla, and splashes from drop impacts on viscous fluids. Physical review letters 122 (1), 014501.
  * Morton et al. (2000) Morton, David, Rudman, Murray & Jong-Leng, Liow 2000 An investigation of the flow regimes resulting from splashing drops. Physics of Fluids 12 (4), 747–763.
  * Mostert et al. (2022) Mostert, Wouter, Popinet, Stéphane & Deike, Luc 2022 High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production. Journal of Fluid Mechanics 942.
  * Murphy et al. (2015) Murphy, David W, Li, Cheng, d’Albignac, Vincent, Morra, David & Katz, Joseph 2015 Splash behaviour and oily marine aerosol production by raindrops impacting oil slicks. Journal of Fluid Mechanics 780, 536–577.
  * Nikolopoulos et al. (2007) Nikolopoulos, Nikos, Theodorakakos, Andreas & Bergeles, George 2007 Three-dimensional numerical investigation of a droplet impinging normally onto a wall film. Journal of computational physics 225 (1), 322–341.
  * Oguz & Prosperetti (1989) Oguz, Hasan N & Prosperetti, Andrea 1989 Surface-tension effects in the contact of liquid surfaces. Journal of Fluid Mechanics 203, 149–171.
  * Oguz & Prosperetti (1990) Oguz, Hasan N & Prosperetti, Andrea 1990 Bubble entrainment by the impact of drops on liquid surfaces. Journal of Fluid Mechanics 219, 143–179.
  * Okawa et al. (2006) Okawa, Tomio, Shiraishi, Takuya & Mori, Toshiaki 2006 Production of secondary drops during the single water drop impact onto a plane water surface. Experiments in fluids 41 (6), 965–974.
  * Okawa et al. (2008) Okawa, Tomio, Shiraishi, Takuya & Mori, Toshiaki 2008 Effect of impingement angle on the outcome of single water drop impact onto a plane water surface. Experiments in Fluids 44 (2), 331–339.
  * Pairetti et al. (2020) Pairetti, Cesar, Damian, Santiago Marquez, Nigro, Norberto, Popinet, Stéphane & Zaleski, Stéphane 2020 Mesh resolution effects on vof simulations of primary atomization. Atomization and Sprays 30 (12), 913–935.
  * Pan et al. (2008) Pan, Kuo-Long, Cheng, Kai-Ren, Chou, Ping-Chung & Wang, Ching-Hua 2008 Collision dynamics of high-speed droplets upon layers of variable thickness. Experiments in fluids 45 (3), 435–446.
  * Philippi et al. (2016) Philippi, Julien, Lagrée, Pierre-Yves & Antkowiak, Arnaud 2016 Drop impact on a solid surface: short-time self-similarity. Journal of Fluid Mechanics 795, 96–135.
  * Popinet (2003) Popinet, Stéphane 2003 Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. Journal of Computational Physics 190 (2), 572–600.
  * Popinet (2009) Popinet, Stéphane 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228 (16), 5838–5866.
  * Popinet (2015) Popinet, Stéphane 2015 A quadtree-adaptive multigrid solver for the serre–green–naghdi equations. Journal of Computational Physics 302, 336–358.
  * Popinet & collaborators (2013–2023) Popinet, Stéphane & collaborators 2013–2023 Basilisk. http://basilisk.fr.
  * Prosperetti et al. (1989) Prosperetti, Andrea, Crum, Lawrence A & Pumphrey, HC 1989 The underwater noise of rain. Journal of Geophysical Research: Oceans 94 (C3), 3255–3259.
  * Prosperetti & Oguz (1993) Prosperetti, Andrea & Oguz, Hasan N 1993 The impact of drops on liquid surfaces and the underwater noise of rain. Annual Review of Fluid Mechanics 25 (1), 577–602.
  * Pumphrey et al. (1989) Pumphrey, Hugh C, Crum, LA & Bjo/Rno/, L 1989 Underwater sound produced by individual drop impacts and rainfall. The Journal of the Acoustical Society of America 85 (4), 1518–1526.
  * Pumphrey & Elmore (1990) Pumphrey, Hugh C & Elmore, Paul A 1990 The entrainment of bubbles by drop impacts. Journal of Fluid Mechanics 220, 539–567.
  * Ray et al. (2015) Ray, Bahni, Biswas, Gautam & Sharma, Ashutosh 2015 Regimes during liquid drop impact on a liquid pool. Journal of Fluid Mechanics 768, 492–523.
  * Reijers et al. (2019) Reijers, Sten A, Liu, Bo, Lohse, Detlef & Gelderblom, Hanneke 2019 Oblique droplet impact onto a deep liquid pool. arXiv preprint arXiv:1903.08978 .
  * Rein (1993) Rein, Martin 1993 Phenomena of liquid drop impact on solid and liquid surfaces. Fluid Dynamics Research 12 (2), 61–93.
  * Rein (1996) Rein, Martin 1996 The transitional regime between coalescing and splashing drops. Journal of Fluid Mechanics 306, 145–165.
  * Resch & Afeti (1991) Resch, François & Afeti, George 1991 Film drop distributions from bubbles bursting in seawater. Journal of Geophysical Research: Oceans 96 (C6), 10681–10688.
  * Resch et al. (1986) Resch, FJ, Darrozes, JS & Afeti, GMhttps 1986 Marine liquid aerosol production from bursting of air bubbles. Journal of Geophysical Research: Oceans 91 (C1), 1019–1029.
  * Riboux & Gordillo (2014) Riboux, Guillaume & Gordillo, José Manuel 2014 Experiments of drops impacting a smooth solid surface: a model of the critical impact speed for drop splashing. Physical review letters 113 (2), 024507.
  * Rieber & Frohn (1999) Rieber, Martin & Frohn, Arnold 1999 A numerical study on the mechanism of splashing. International Journal of Heat and Fluid Flow 20 (5), 455–461.
  * Rioboo et al. (2002) Rioboo, R, Marengo, M & Tropea, C 2002 Time evolution of liquid drop impact onto solid, dry surfaces. Experiments in fluids 33 (1), 112–124.
  * Rodriguez & Mesler (1985) Rodriguez, Francisco & Mesler, Russell 1985 Some drops don’t splash. Journal of colloid and interface science 106 (2), 347–352.
  * Rudman (1998) Rudman, Murray 1998 A volume-tracking method for incompressible multifluid flows with large density variations. International Journal for numerical methods in fluids 28 (2), 357–378.
  * Samet (1990) Samet, Hanan 1990 The design and analysis of spatial data structures, , vol. 85. Addison-Wesley Reading, MA.
  * Sanjay et al. (2022) Sanjay, Vatsal, Lakshman, Srinath, Chantelot, Pierre, Snoeijer, Jacco H & Lohse, Detlef 2022 Drop impact on viscous liquid films. arXiv preprint arXiv:2206.06298 .
  * Schotland (1960) Schotland, RM 1960 Experimental results relating to the coalescence of water drops with water surfaces. Discussions of the Faraday Society 30, 72–77.
  * Shin et al. (2017) Shin, Seungwon, Chergui, Jalel & Juric, Damir 2017 A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows. Journal of Mechanical Science and Technology 31 (4), 1739–1751.
  * Sochan et al. (2018) Sochan, Agata, Beczek, Michał, Mazur, Rafał, Ryżak, Magdalena & Bieganowski, Andrzej 2018 The shape and dynamics of the generation of the splash forms in single-phase systems after drop hitting. Physics of Fluids 30 (2), 027103\.
  * Thomson & Newall (1886) Thomson, Joseph John & Newall, Hugh Frank 1886 V. on the formation of vortex rings by drops falling into liquids, and some allied phenomena. Proceedings of the royal society of London 39 (239-241), 417–436.
  * Thoraval et al. (2016) Thoraval, Marie-Jean, Li, Yangfan & Thoroddsen, Sigurdur T 2016 Vortex-ring-induced large bubble entrainment during drop impact. Physical Review E 93 (3), 033128.
  * Thoraval et al. (2013) Thoraval, M-J, Takehara, K, Etoh, TG & Thoroddsen, Sigurdur T 2013 Drop impact entrapment of bubble rings. Journal of Fluid Mechanics 724, 234–258.
  * Thoraval et al. (2012) Thoraval, Marie-Jean, Takehara, Kohsei, Etoh, Takeharu Goji, Popinet, Stéphane, Ray, Pascal, Josserand, Christophe, Zaleski, Stéphane & Thoroddsen, Sigurdur T 2012 von kármán vortex street within an impacting drop. Physical review letters 108 (26), 264506.
  * Thoroddsen (2002) Thoroddsen, ST 2002 The ejecta sheet generated by the impact of a drop. Journal of Fluid Mechanics 451, 373.
  * Thoroddsen et al. (2003) Thoroddsen, ST, Etoh, TG & Takehara, K 2003 Air entrapment under an impacting drop. Journal of Fluid Mechanics 478, 125–134.
  * Thoroddsen et al. (2009) Thoroddsen, Sigurdur T, Takehara, K, Etoh, TG & Ohl, C-D 2009 Spray and microjets produced by focusing a laser pulse into a hemispherical drop. Physics of Fluids 21 (11), 112101\.
  * Thoroddsen et al. (2011) Thoroddsen, Sigurdur T, Thoraval, M-J, Takehara, K & Etoh, TG 2011 Droplet splashing by a slingshot mechanism. Physical review letters 106 (3), 034501.
  * Thoroddsen et al. (2012) Thoroddsen, Sigurdur T, Thoraval, M-J, Takehara, K & Etoh, TG 2012 Micro-bubble morphologies following drop impacts onto a pool surface. Journal of Fluid Mechanics 708, 469–479.
  * Tran et al. (2013) Tran, Tuan, de Maleprade, Hélène, Sun, Chao & Lohse, Detlef 2013 Air entrainment during impact of droplets on liquid surfaces. Journal of Fluid Mechanics 726.
  * Veron (2015) Veron, Fabrice 2015 Ocean spray. Annu. Rev. Fluid Mech 47 (1), 507–538.