1. Introduction
The study of fluid drop impacts on a liquid surface is fundamental to engineering applications such as spray cooling and inkjet printing (Liang & Mudawar Reference Li, Schmitt, Ummenhofer and Karnauskas2016). Furthermore, similar processes in nature, particularly rainfall on the ocean, play an important role in influencing ocean surface salinity and temperature (Li et al. Reference Li, Jia, Liu and Su2016b ; Balaguru et al. Reference Balaguru, Foltz, Leung and Hagos2022; Shackelford et al. Reference Schiller and Naumann2022). Ocean raindrops can also bias rainfall measurements, as instruments may overestimate rainfall intensity by detecting secondary droplets and sprays (Okachi et al. Reference Oguz and Prosperetti2020), contribute to oil spill dispersion (Wu et al. Reference Wu, Hannah, Thupaki, Mo and Law2017), and generate acoustic signals that can be used to monitor precipitation patterns (Prosperetti & Oguz Reference Poureslami, Faghiri and Shafii1993).
Worthington (Reference Worthington1883) pioneered the visualisation of drop–surface interactions, laying the foundation for the extensive research that followed. Many configurations have since been studied, including impacts on both solid and liquid surfaces, as well as gas film-mediated dynamics (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023; Sprittles Reference Sochan, Beczek, Mazur, Ryżak and Bieganowski2024). For liquid surfaces, the dynamics differs markedly between shallow and deep pools: shallow pools exhibit confined cavity formation and suppressed splashing due to proximity to the solid base (Li et al. Reference Leng2019; Tang et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2019; Wu & Saha Reference Wu and Saha2022), whereas deep pools permit larger cavities, higher jets, and more vigorous splashing (Sykes et al. Reference Sprittles2023). Murphy et al. (Reference Morton, Rudman and Liow2015) categorised previous studies on drop impact on a deep pool into five regimes based on their Froude number
$\textit{Fr} = U_0^2/(g d)$
and Weber number
$\textit{We} = \rho _l U_0^2 d/\sigma$
, where
$U_0$
is the drop impact velocity,
$g$
is gravity,
$d$
is the drop diameter,
$\rho _l$
is the liquid density, and
$\sigma$
is the surface tension. The first regime, coalescence and vortex ring (C&VR,
$40\lt \textit{Fr}\lt 100$
,
$\textit{We}\lt 100$
), involves drops that fall with negligible momentum and coalesce with the liquid, where high surface tension results in the formation of a downward-moving vortex ring after impact (Hsiao et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel1988; Deng et al. Reference Deng, Anilkumar and Wang2007; Dong & Angeli Reference Dong and Angeli2023; Anirudh et al. Reference Anirudh, Behera and Sahu2024; Behera et al. Reference Behera, Deka, Sahu and Biswas2026). The second regime, swell and thin jet (S&TJ,
$100\lt \textit{Fr}\lt 500$
,
$100\lt \textit{We}\lt 300$
), is characterised by the formation of an impact cavity and outward swell; the cavity may retract and form a rising liquid jet that disintegrates into secondary droplets (Manzello & Yang Reference Liu, Zhang, Zheng and Duncan2002; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2015). The regular entrainment regime is a sub-regime of S&TJ, where capillary waves travel down the cavity surface and pinch off an air bubble (Elmore et al. Reference Elmore, Chahine and Oguz2001; Leng Reference Krechetnikov and Homsy2001). At higher energies, the crown and thick jet regime (
$100\lt \textit{Fr}\lt 1000$
,
$300\lt \textit{We}\lt 2000$
) features increased cavity depth, while the outward swell rises to form a crown, from which ligaments form and disintegrate into secondary droplets. The rebounding cavity also forms a thick rising liquid jet, and the initial impact generates a cloud of small drops (Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2015). The most energetic regime, bubble canopy (BC,
$300\lt \textit{Fr}\lt 7000$
,
$2000\lt \textit{We}\lt 20\,000$
), produces large cavities and crowns; the crown rim rises and expands radially before converging to form a canopy and a liquid jet that travels in both vertical directions. The downward-travelling jet can penetrate the cavity bottom and entrain bubbles into the liquid pool (Murphy et al. Reference Morton, Rudman and Liow2015; Sochan et al. Reference Shackelford, DeMott, van Leeuwen, Thompson and Hagos2018; Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). Previous studies on multi-drop impact on pools have largely focused on the C&VR and S&TJ regimes (Santini et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2017; Guilizzoni et al. Reference Guilizzoni, Santini and Fest-Santini2019; Guilizzoni & Frontera Reference Guilizzoni and Frontera2022; Poureslami et al. Reference Popinet2023; Wang et al. Reference Wang, Li, Chong and Zhou2024). In comparison, very little attention has been given to multi-drop interactions in the BC regime, where the impact morphology becomes highly complex, and large numbers of secondary droplets are produced within a very short period. Since raindrops fall in the BC regime (Gunn & Kinzer Reference Gunn and Kinzer1949; Murphy et al. Reference Morton, Rudman and Liow2015), it is important to understand how neighbouring raindrop impacts interact, and how these interactions influence the generation, size distribution and dynamics of secondary droplets.
Secondary droplets generated during drop impacts originate from several mechanisms. Prompt splash, which results from instabilities in the initial ejecta sheet, produces very fine droplets within a very short time before the drop fully merges with the pool (Deegan et al. Reference Deegan, Brunet and Eggers2007; Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012; Marcotte et al. Reference Manzello and Yang2019). Crown rim break-up occurs later as rim ligaments fragment into medium and large droplets (Villermaux et al. Reference Villermaux, Marmottant and Duplat2004; Eggers & Villermaux Reference Eggers and Villermaux2008; Murphy et al. Reference Morton, Rudman and Liow2015). In multi-drop impacts, a distinguishing feature is the central liquid film, a vertically rising sheet that forms between adjacent droplets upon impact. Previous studies have shown that this film can reduce the splashing threshold and provide an additional pathway for secondary droplet formation (Liang et al. Reference Liang, Zhang, Chen, Yu and Shen2018; Fest-Santini et al. Reference Fest-Santini, Steigerwald, Santini, Cossali and Weigand2021; Zhou et al. Reference Zhou, Jin, Yang and Yu2024). Studies of drop impact on solid substrates have also examined secondary droplet dynamics and interaction effects, including analytical models for the diameters and velocities of ejected droplets (Riboux & Gordillo Reference Rayleigh2015), data-driven analyses of splashing thresholds (Pierzyna et al. Reference Pairetti, Damian, Nigro, Popinet and Zaleski2021), quantification of atomisation in prompt splash and crown rim break-up (Burzynski et al. Reference Burzynski, Roisman and Bansmer2020), and interactions involving multiple impacting drops and non-isolated impacts on substrates (Roisman et al. Reference Roisman, Horvat and Tropea2002; Wang & Bourouiba Reference Wang and Bourouiba2018). In the context of secondary droplet formation during high-speed impacts representative of rainfall, Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023) examined the droplet population up to
$4$
ms after impact, and emphasised the generation of very small droplets associated with the prompt splash. Liu et al. (Reference Liang, Zhang, Yu, Chen and Shen2024) recreated rainfall conditions in the laboratory, and measured the secondary droplet number density at various heights above the surface. However, the highly dynamic nature of high-speed impacts has made it difficult to quantitatively determine the mechanisms responsible for secondary droplet formation, leaving a substantial gap in current understanding (Veron Reference Veron2015).
Simulation-based studies on drop–pool interactions in the literature focus their efforts on resolving the interfacial dynamics between different phases. Oguz & Prosperetti (Reference Murphy, Li, d’Albignac, Morra and Katz1990) used a boundary integral model to simulate drop impact, which included surface tension effects at the interface through curvature, with stabilisation to keep the surface smooth. However, this method assumes irrotational flow, which fails to predict vortex rings in the pool. Morton et al. (Reference Mathai, Lohse and Sun2000) addressed this by utilising the continuum surface force model proposed by Brackbill et al. (Reference Brackbill, Kothe and Zemach1992), where surface tension is added as a body force in the momentum equation, with the force direction given by the interface normal, and its strength determined by the interface curvature. In terms of simulation construction, two-dimensional (2-D) simulations on a planar domain (Coppola et al. Reference Coppola, Rocco and de Luca2011; Kharmiani et al. Reference Hsiao, Lichter and Quintero2016) significantly reduce computational costs, but provide only partial insight into the drop impact process in the radial direction. Additionally, they physically represent an infinite cylinder instead of a spherical drop. To address this, 2-D axisymmetric simulations are more commonly employed while still keeping computational costs low (Fudge et al. Reference Fudge, Cimpeanu and Castrejón-Pita2021, Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023; Sanjay et al. Reference Roisman, Prunet-Foch, Tropea and Vignes-Adler2023; Wang et al. Reference Wang, Li, Chong and Zhou2024, Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2025). To capture non-axisymmetric features, three-dimensional (3-D) simulations have been performed in a computational domain that is restricted to half or a quarter of the drop geometry, where the symmetric boundary condition is applied (Brambilla & Guardone Reference Brambilla and Guardone2015; Guo & Lian Reference Guo and Lian2017; Guilizzoni et al. Reference Guilizzoni, Santini and Fest-Santini2019; Guilizzoni & Frontera Reference Guilizzoni and Frontera2022). Additionally, 3-D simulations often employ adaptive mesh refinement (AMR) to reduce the overall mesh count while preserving sufficient resolution, as demonstrated by Cimpeanu & Papageorgiou (Reference Cimpeanu and Papageorgiou2018), Liang et al. (Reference Liang, Zhang, Chen, Yu and Shen2018, Reference Liang and Mudawar2019) and Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). Since axisymmetric 2-D simulations or 3-D simulations with symmetric domains may exclude transverse instabilities in the crown and the resulting secondary droplet shedding (Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2025), which are crucial to the overall impact dynamics, a fully 3-D simulation is more suitable for investigating the impact morphology and secondary droplet formation.
The present study addresses these gaps by constructing a fully 3-D numerical simulation of high-speed raindrops impacting a liquid pool to represent rainfall on the ocean. A dimensional analysis is performed to quantify the number density of secondary droplets as a function of the droplet radius, with the aim of providing insight into the key factors governing the droplet generation process. The impact morphology of raindrops is examined, and the stages of their interaction identified and compared across multiple inter-raindrop spacings. The spatial and temporal distributions of secondary droplets are then analysed to reveal how the fate of these droplets is shaped by raindrop interactions. The remainder of this paper is organised as follows. The numerical solver, simulation set-up and validation are described in § 2. In § 3, we present the results and discussion, including early-time secondary droplet size distributions, impact morphology, and the influence of raindrops on secondary droplet evolution. Conclusions are given in § 4.
2. Methodology
2.1. Numerical solver
The current study carries out direct numerical simulations of the incompressible air–water system using the open-source solver Basilisk (Popinet Reference Pierzyna, Burzynski, Bansmer and Semaan2018). In Basilisk, the volume of fluid method is used to track the evolution of the interface by following the split advection method (Weymouth & Yue Reference Weymouth and Yue2010) to solve the advection equation
\begin{equation} \frac {\partial f}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u} f) = 0, \end{equation}
where
$f$
represents the volume fraction of water in each cell. The Navier–Stokes equations with surface tension are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0, \nonumber \\ \rho \frac {\partial \boldsymbol{u}}{\partial t} + \rho (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u} &= -\boldsymbol{\nabla }\!p + \boldsymbol{\nabla }\boldsymbol{\cdot } \big[\mu (\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^{\rm T})\big] + \rho \boldsymbol{g} + \sigma \kappa \delta _s \boldsymbol{n}, \end{align}
where
$\boldsymbol{u}$
is the velocity vector,
$\rho =f\rho _{w}+(1-f)\rho _{a}$
is the density,
$\mu =f\mu _{w}+ (1-f)\mu _{a}$
is the dynamic viscosity (subscripts
$a$
and
$w$
represents air and water, respectively),
$p$
is the pressure,
$\boldsymbol{g}$
is gravity,
$\sigma$
is the surface tension coefficient,
$\kappa$
is the curvature of the air–water interface,
$\delta _s$
is a surface Dirac function that is non-zero on the interface and zero everywhere else, and
$\boldsymbol{n}$
is the unit vector normal to the interface. For spatial discretisation, Basilisk employs an octree-structured grid, in which each cubic cell is divided into eight smaller cells at each refinement level. The current study also utilises the AMR function in Basilisk (van Hooft et al. Reference Thoroddsen2018), which dynamically refines or coarsens the grid based on the discretisation error in spatial fields, i.e.
$f$
and
$\boldsymbol{u}$
. This error is compared with a predefined threshold to determine whether a cell should be refined, coarsened, or left unchanged. More details of the numerical scheme used in Basilisk can be found in Popinet (Reference Pierzyna, Burzynski, Bansmer and Semaan2018).
2.2. Simulation set-up
In the present study, a series of raindrop impact simulations with different configurations is performed. The single-raindrop simulations with different surface tension values and diameters are labelled as SR, SR-
$0.5\sigma _0$
, SR-
$0.25\sigma _0$
, SR-
$2\sigma _0$
, SR-
$0.25d_0$
and SR-
$0.5d_0$
, with
$\sigma _0$
and
$d_0$
denoting the surface tension and the effective raindrop diameter of case SR, respectively. Case SR serves as the reference case for comparison with all other simulations, while the remaining single-raindrop cases are performed to investigate the effect of surface tension and raindrop size on early secondary droplet formation. In addition, three two-raindrop simulations are performed with different centre-to-centre distances in the
$x$
-direction:
$D=2d_h, 3d_h, 4d_h$
, denoted as D2, D3, D4, respectively. They are designed to investigate how raindrop interaction influences impact morphology and secondary droplet dynamics. Because the single-raindrop cases (except for SR) focus on early secondary droplet break-up, they are only simulated until 5.7 ms in physical time, whereas in cases SR, D2 and D3, the simulations are run to 45.6 ms, and D4 is extended to 47.6 ms to fully capture the interface evolution. Figure 1 shows a schematic of the layout of the two-raindrop simulation domain. More details on the parameters used in each simulation case can be found in table 1.
Physical parameters used in the simulations. Here,
$N_{\textit{rain}}$
is the number of raindrops,
$U_0$
is the initial raindrop velocity,
$d_h$
and
$d_v$
are the major and minor axes of the ellipsoidal raindrop,
$\sigma /\sigma _0$
is the surface tension normalised by the reference value
$\sigma _0$
of case SR, and
$D$
is the inter-raindrop distance.
Table 1 Long description
The table presents a comparison of different raindrop impact simulations with various configurations. It includes 12 rows and 10 columns. The columns are labeled as Case, N rain, L max, U 0 (m s^-1), d h (mm), d v (mm), σ/σ 0, D/d h, and Simulation duration (ms). The rows list different simulation cases with their respective parameters. For example, Row 1: Case SR, N rain 1, L max 12, U 0 7.2, d h 4.3, d v 3.8, σ/σ 0 1, D/d h 1, Simulation duration 45.6 ms. Row 2: Case SR-L10, N rain 1, L max 10, U 0 7.2, d h 4.3, d v 3.8, σ/σ 0 1, D/d h not applicable, Simulation duration 45.6 ms. The table provides detailed parameters for each simulation case, including the number of raindrops, initial raindrop velocity, major and minor axes of the ellipsoidal raindrop, normalized surface tension, and inter-raindrop distance.
Layout of the computational domain with a cross-sectional view in the
$x$
–
$y$
plane illustrating the initial configuration of the raindrops and free surface in a two-raindrop case.
Figure 1 Long description
Panel A: A 3D diagram showing the computational domain with air above and water below. Two raindrops are positioned above the water surface, ready to impact. Panel B: A cross-sectional view in the x-y plane detailing the initial configuration. The diagram includes labels for the free surface, air, water, and the initial mesh. The raindrops are shown with their diameters and the distance between them. The velocity of the raindrops and the depth of the water are also indicated.
The physical parameters of the simulations are similar to those in the laboratory experiment of Murphy et al. (Reference Morton, Rudman and Liow2015) and the numerical study of Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). The impacting raindrops are ellipsoids, with major axis length
$d_h=4.3\text{ mm}$
and minor axis length
$d_v=3.8\text{ mm}$
. The effective diameters are defined as
$d_0=(d_vd_h^2)^{1/3}=4.1\text{ mm}$
. Raindrops of this size fall within the less-studied BC regime, which characterises the most energetic air–water impact dynamics, and occur during intense rainfall (Best Reference Best1950). The initial raindrop velocity before impact is set to
$U_0 =7.2 \rm\,{ m\,s}^{-1}$
, which is approximately 82 % of the terminal velocity of raindrops of this size falling in quiescent air based on the field observations from Gunn & Kinzer (Reference Gunn and Kinzer1949), which show that the velocity of the raindrops varies in the range 0.2–9 m s−1, and their diameters typically range from approximately 0.1 mm to 6 mm. The initial gap between the raindrop and the free surface is set to
$\delta =0.1d_0$
.
As illustrated in figure 1, the computational domain length is selected as
$L=16d_0$
to ensure a balance between computational efficiency and minimisation of boundary effects. The pool depth is set to
$H=8d_0$
, providing sufficient space above the free surface to form secondary droplets from splashing. An outflow boundary condition is applied at the top of the domain, while all other boundaries use a symmetry boundary condition. The density and viscosity ratios between water and air are
$\rho _w/\rho _a=783$
and
$\mu _w/\mu _a=56$
, respectively. The following dimensionless numbers can be obtained: Reynolds number
$\textit{Re} = \rho _{w} U_0 d_0/\mu _{w} = 30\,069$
, Weber number
$\textit{We} = \rho _w U_0^2 d_0/\sigma _0 = 2964$
, and Froude number
$\textit{Fr} = U_0^2/(g d_0) = 1290$
. The computational grid is initialised based on the distance to the droplet centres: the mesh is refined up to level 12 (corresponding to a smallest mesh size
$\varDelta = 15.6\,{\unicode{x03BC}} \text{m}$
) within the raindrops, and progressively coarsened with increasing distance. After initialisation, the AMR algorithm adaptively refines the grid, with maximum refinement level
$L_{\textit{max}}=12$
. As an illustration of the computational cost of our simulations, case D2 consumes approximately
$7\times 10^4$
core-hours to simulate 45.6 ms of the impact process.
2.3. Validation
For validation, simulation data from the single-raindrop case SR are used to calculate four quantities representative of the impact morphology: the cavity radius
$r_c$
, defined as the radius of the circular interface at
$y=0$
, the cavity depth
$D_c$
, defined as the distance between the deepest point of the cavity interface and the free surface at
$y=0$
, the crown rim height
$H_{\textit{rim}}$
, defined as the distance between the base of the ligaments and the free surface at
$y=0$
, and the crown rim radius
$r_{\textit{rim}}$
, defined as the radial distance from the impact centre to the crown rim at the identified rim height. The full procedure for extracting
$r_c$
,
$D_c$
,
$H_{\textit{rim}}$
and
$r_{\textit{rim}}$
is detailed in Appendix A. To assess grid convergence and grid independence, the time histories of
$r_c$
,
$D_c$
,
$H_{\textit{rim}}$
and
$r_{\textit{rim}}$
for the SR case are compared using
$L_{\textit{max}}=10$
, 11 and 12, as shown in figure 2. While the results for
$L_{\textit{max}}=10$
and 11 deviate from those reported by Murphy et al. (Reference Morton, Rudman and Liow2015) and Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023), the results for
$L_{\textit{max}}=12$
show good agreement with these reference studies, suggesting that
$L_{\textit{max}}=12$
is sufficient for capturing the impact morphology in the present study.
In figure 3, we compare the secondary droplet size distribution averaged over
$0.57$
–
$4\,\mathrm{ms}$
with results obtained over
$0.2$
–
$4\,\mathrm{ms}$
by Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023), who used a higher resolution. A notable improvement in our distribution is observed when
$L_{\textit{max}}$
is increased from
$11$
to
$12$
. The remaining difference between our results and theirs is primarily attributed to temporal variability in the droplet production: for instance, their results have shown that the droplet number can vary by up to one order of magnitude between two time windows
$0.2$
–
$2\,\mathrm{ms}$
and
$3$
–
$4\,\mathrm{ms}$
. Overall, this comparison suggests that the resolution used in the present study is sufficient to capture the dominant secondary droplet population.
Time evolution of (a) cavity radius, (b) cavity depth, (c) crown rim height, and (d) crown rim radius in case SR for
$L_{\textit{max}}=10$
, 11 and 12. Also plotted are reference results from Murphy et al. (Reference Morton, Rudman and Liow2015) and Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023).
Time-averaged secondary droplet size distributions extracted from case SR over 0.57–4 ms for
$L_{\textit{max}}=11$
and
$12$
, compared with the results of Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023) averaged over 0.2–4 ms.
3. Results
3.1. Early generation of secondary droplets
In this subsection, we examine the early generation of secondary droplets from raindrop impacts. The focus is on the early stage of raindrop–pool interaction within
$5.7\ \mathrm{ms}$
after release, when most secondary droplets are produced. Following the discussion of mesh resolution on atomisation and droplet generation from Pairetti et al. (Reference Okachi, Yamada, Baba and Kubo2020), the droplet statistics in the present study only include droplets with diameters larger than
$64\,{\unicode{x03BC}} \text{m}$
, approximately four times the size of the smallest cell. Secondary droplets smaller than this value are considered not well resolved for our mesh set-up, and are not discussed here.
3.1.1. Possible mechanism for secondary droplet production
Historically, the Rayleigh–Plateau instability (Rayleigh Reference Rayleigh1892; Bremond & Villermaux Reference Bremond and Villermaux2006), the Richtmyer–Meshkov instability (Richtmyer Reference Riboux and Gordillo1960; Gueyffier & Zaleski Reference Gueyffier and Zaleski1998), and a combined Richtmyer–Meshkov and Rayleigh–Taylor mechanism (Rayleigh Reference Prosperetti and Oguz1882; Krechetnikov & Homsy Reference Kolmogorov2009) have been proposed as possible candidates for governing droplet generation from the crown rim. Since previous studies have attributed secondary droplet production during drop–pool impact to Rayleigh–Plateau instability of the crown rim (Roisman et al. Reference Richtmyer2006; Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita, Matar and Castrejón-Pita2023; Anirudh et al. Reference Anirudh, Behera and Sahu2026), we first examine whether this mechanism can explain the secondary droplet generation observed here. Following Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010), we compare the measured crown rim corrugation spectrum with the Rayleigh–Plateau prediction based on the instantaneous crown radius and rim thickness, as described in Appendix B. Our calculations indicate that crown break-up in the present high-
$\textit{Re}$
and high-
$\textit{We}$
regime is not governed by direct Rayleigh–Plateau growth on the crown rim.
Using the impacting raindrop diameter and velocity in the present simulations gives
$\textit{Re}_p=U_0d_0/\nu _a\approx 2000$
, placing the airflow around the falling raindrop in the unsteady and turbulent wake regime (Mathai et al. Reference Marcotte, Michon, Séon and Josserand2020). Thus the early impact and crown development occur in a strongly inertial gas-flow environment, in which the break-up process is likely affected by unsteady airflow. In the fragmentation framework reviewed by Villermaux (Reference Villermaux2007), the log-normal distribution is associated with a sequential cascade break-up mechanism (Kolmogorov Reference Kharmiani, Passandideh-Fard and Niazmand1941), whereas gamma-type distributions are linked to ligament-mediated break-up and the coalescence and rearrangement of corrugations along ligaments (Bremond & Villermaux Reference Bremond and Villermaux2006). Following these ideas, we have examined whether the normalised droplet radius distributions could be represented by log-normal or gamma-type forms. However, neither provides a robust description across the cases and times considered (not shown). In fact, the difference between these two distributions in quantifying droplets generated during fragmentation is relatively small (Villermaux Reference Villermaux2007). Therefore, distinguishing which of these two mechanisms governs the present problem will require additional simulations to obtain droplet statistics with reduced scatter in future studies.
3.1.2. Scaling law derivation and validation
We use the Kolmogorov–Hinze framework (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955) to guide the scaling analysis of the secondary droplet size distribution. This approach provides a dimensional argument for how droplet production may be controlled by the competition between inertial fluctuations and surface tension. In the present analysis, it is used to describe the surface-tension-controlled part of the spectrum. We do not attempt to derive a separate inertia-dominated scaling for larger droplets because the sample size is limited and the scatter is larger. One notable application of this framework is found in studies of air bubbles in breaking waves, where different scaling laws were identified for bubbles larger (Garrett et al. Reference Garrett, Li and Farmer2000) and smaller (Deane & Stokes Reference Deane and Stokes2002) than the Hinze scale. Inspired by these seminal works, we derive a power scaling law for the size spectrum of secondary droplets
$N_d$
, i.e. the number of secondary droplets per unit droplet radius increment per unit volume of air. It is assumed that
$N_d$
is a multiplicative function of the rate of water being ejected into the air
$Q$
– i.e. the volume of water ejected per volume of air per second – surface tension
$\sigma$
, water density
$\rho _w$
, and droplet radius
$r_s$
. After equalising the dimensions of these parameters, the following relationship between
$N_d(r_s)$
and the droplet radius
$r_s$
can be obtained:
\begin{equation} N_d(r_s)\propto Q(\sigma /\rho _w)^{-1/2}r_s^{-5/2}. \end{equation}
It is worth noting that while (3.1) is similar to the relationship derived by Deane & Stokes (Reference Deane and Stokes2002), we do not introduce the raindrop impact velocity as an additional independent parameter in the present scaling. Here,
$Q$
is treated as the impact-driven water ejection rate, and can therefore vary with the impact conditions, including the impact velocity. Thus the present scaling should not be interpreted as assuming that velocity effects are negligible; rather, the effect of velocity is incorporated through its influence on
$Q$
.
Time-averaged number density spectra as functions of the droplet radius, normalised by the number of raindrops
$N_{\textit{rain}}$
for cases SR, D2, D3 and D4.
To validate (3.1), the time-averaged secondary droplet number density spectra are first computed for cases SR, D2, D3 and D4 over the interval
$t = 0.57$
ms to
$t = 5.7$
ms. Figure 4 shows the result normalised by the number of raindrops
$N_{\textit{rain}}$
. The collapse of the data among different cases confirms a negligible effect associated with the raindrop interaction. More importantly, all four cases follow the predicted
$N_d(r_s) \propto r_s^{-5/2}$
scaling well.
(a) Time-averaged number density spectra of secondary droplets, and (b) the corresponding rescaled spectra as functions of droplet radius for single-raindrop cases with different surface tensions.
Same as figure 5 but for single-raindrop cases with different raindrop diameters.
A similar strategy is adopted for investigating the impact of surface tension on the secondary droplet distributions. Their time-averaged spectra are shown in figure 5(a), demonstrating the impact caused by varying the surface tension. Figure 5(b) shows the rescaled spectra
$N_d(r_s)(\sigma _0/\sigma )^{-1/2}$
, which collapse well in the range
$0.05\lt r_s \lt 0.1\ \mathrm{mm}$
for the cases with different surface tension, thus confirming the scaling, i.e.
$N_d(r_s) \propto \sigma ^{-1/2}$
. Finally, we examine how the raindrop size affects the secondary droplet distribution. Figure 6(a) shows the time-averaged spectra of the single-raindrop simulations with varying raindrop diameters, including cases SR, SR-
$0.5d_0$
and SR-
$0.25d_0$
(see table 1). All cases follow the
$N_d(r_s) \propto r_s^{-5/2}$
scaling closely, indicating that the underlying break-up mechanism remains similar across different raindrop sizes.
For different raindrop diameters, the secondary droplet number density continues to follow the same
$r_s^{-5/2}$
scaling. The main effect of increasing
$d$
is instead an increase in the magnitude of
$N_d$
. As shown in figure 6(b), the spectra collapse well when the number densities are multiplied by
$(d_0/d)^2$
, suggesting a dependence
$N_d\propto d^2$
. This scaling can be interpreted from the crown volume: for a thin crown sheet, the volume is expected to scale approximately with the raindrop impact footprint, which is proportional to
$d^2$
. The resulting increase in crown volume provides more liquid available for break-up into secondary droplets, consistent with the collapse of the crown volume time histories in figure 7.
Time histories of the non-dimensionalised crown volume rescaled by
$(d_0/d)^2$
for single-raindrop cases with different drop diameters.
We briefly discuss the scatter observed in the simulation data, and its possible causes. For example, the smallest secondary droplets approach the under-resolved size limit of the mesh, which can lead to over-counting, as noted by Pairetti et al. (Reference Okachi, Yamada, Baba and Kubo2020). As expected, this numerical issue becomes more significant as the surface tension decreases (comparing
$N_d$
for
$r_s \lt 0.04\space \mathrm{mm}$
in figures 5
a,b). On the other hand, the scatter also tends to increase for droplets larger than approximately
$0.1$
mm because fewer larger droplets are generated. For instance, there are approximately
$O(10)$
secondary droplets with a radius of
$0.2$
mm in the computational domain. Because of the scatter, we do not attempt to derive a scaling law for droplets larger than the Hinze scale, for which break-up induced by inertial fluctuations is expected to dominate over surface tension.
3.2. Impact morphology
In this subsection, the morphology of the impact is analysed. As the single-raindrop case has been extensively studied, the focus here is on the two-raindrop cases. For all cases,
$t = 0$
is defined as the start of the simulation, when the raindrops are initially positioned above the free surface (figure 1). Figures 8, 9 and 10 illustrate several stages of the impact process for case D2, D3 and D4. First, upon impact, the raindrops create flat-bottomed cavities resembling those seen in single-raindrop impacts. These cavities gradually transition into hemispherical shapes while producing cylindrical upward liquid sheets, or crowns (Thoroddsen Reference Tang, Saha, Law and Sun2002; Deegan et al. Reference Deegan, Brunet and Eggers2007; Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012). Second, the crowns merge, forming a central liquid film, a feature observed in previous multi-drop impact studies (Li et al. Reference Li, Zhang and Liu2016a
; Liang et al. Reference Liang, Zhang, Chen, Yu and Shen2018; Fest-Santini et al. Reference Fest-Santini, Steigerwald, Santini, Cossali and Weigand2021; Poureslami et al. Reference Popinet2023; Zhou et al. Reference Zhou, Jin, Yang and Yu2024). Third, the subsurface cavities interact and eventually merge into a single large cavity after the central film fully ruptures and detaches from the cavity wall. Finally, the release of surface tension following the detachment of the central film from the cavity bottom leads to a recoil that generates a localised downward bulge at the cavity bottom, further deepening the cavity. Notably, in high-speed single-raindrop impacts, it is common for the crown to converge radially inwards and pinch off, creating a BC (Worthington Reference Worthington1883; Engel Reference Engel1966; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Morton, Rudman and Liow2015). In the two-raindrop cases (D2, D3, D4), although the crowns also curve inwards and partially merge, they do not fully enclose at the centre. As a result, the characteristic jetting behaviour observed after BC formation in single-raindrop cases is absent here.
(a–d) Side views and (e–h) isometric views of the impact morphology in case D2. In (e–h), only part of the domain is shown for clarity. The time instants in (a–d) are 1.14, 6.27, 17.67 and 27.65 ms, respectively, whereas those in (e–h) are 8.83, 11.68, 17.67 and 27.65 ms, respectively.
Same as figure 8, but for case D3. The time instants in (a–d) are 1.14, 7.98, 29.93 and 41.33 ms, respectively, whereas those in (e–h) are 25.08, 27.36, 29.93 and 41.33 ms, respectively.
Comparing the side views shown in figures 8, 9 and 10, it is clear that increasing the inter-raindrop distance progressively delays the merging of crowns and cavities, yielding a weaker downward bulge and shallower penetration depth. This trend is closely linked to the central film behaviour. With larger spacing, the film disintegration (depicted on the bottom rows of figures 8, 9 and 10) occurs later and with less intensity, which reduces the recoil upon detachment from the cavity base, and suppresses bulge formation.
Same as figure 8, but for case D4. The time instants in (a–d) are 1.14, 11.4, 45.03 and 47.6 ms, respectively, whereas those in (e–h) are 22.8, 33.92, 45.03 and 47.6 ms, respectively.
To examine the temporal evolution of cavity depth, we first outline the analytical solution for a single-raindrop impact formulated by Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010). In this model, based on potential flow theory, the velocity potential of the downward-expanding cavity is treated as a superposition of flow past a sphere and a radially expanding flow. The pressure distribution at the cavity surface is determined using the Bernoulli equation, and a dynamic boundary condition is applied to account for capillary forces and gravity. These lead to a system of ordinary differential equations governing the evolution of the cavity geometry:
\begin{eqnarray} \ddot {\alpha } & = & -\frac {3 \dot \alpha ^2}{2 \alpha } - \frac {2}{\alpha ^2\, \textit{We}} - \frac {1}{\textit{Fr}}\frac {\zeta }{\alpha } + \frac {7 \dot \zeta ^2}{4 \alpha }, \end{eqnarray}
\begin{eqnarray} \ddot {\zeta } & = & -3\frac { \dot \alpha \dot {\zeta }}{\alpha } - \frac {9}{2}\frac { \dot \zeta ^2}{ \alpha } - \frac {2}{\textit{Fr}}, \end{eqnarray}
where
$\alpha$
and
$\zeta$
are the cavity radius and the axial position of the cavity centre, normalised by the initial droplet diameter. The initial conditions for this model are derived from the pressure distribution at the impact axis and the free surface of the cavity, resulting in
\begin{equation} \dot {\alpha } \approx 0.17, \quad \alpha \approx \alpha _0 + 0.17 \tau , \quad \dot {\zeta } \approx 0.27, \quad \zeta \approx -\alpha _0 + 0.27 \tau , \end{equation}
where
$\tau =U_0t/d_0$
is the non-dimensional time, and
$\alpha _0$
is a constant associated with the cavity radius
$r_c$
. To determine
$\alpha _0$
, one can use its connection to the dimensionless cavity diameter
$\varOmega =2r_c/d_0$
:
\begin{equation} \varOmega = 2 \sqrt {\alpha ^2 - \zeta ^2} \approx 2 \sqrt {(\alpha _0 + 0.17 \tau )^2 - (0.27 \tau - \alpha _0)^2}. \end{equation}
By fitting the cavity radius data gathered from our SR simulation (refer to figure 2
a),
$\alpha _0\approx 0.81$
is obtained. Note that in the derivation by Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010), the initial condition is determined to be at
$\tau =2$
because it marks the transition from complex interaction between the droplet and pool surface to single-fluid inertial flow, which is more suitable for potential flow analysis. Therefore, the initial conditions for (3.2) and (3.3) are
$ \alpha (2) = 1.15$
,
$ \dot {\alpha }(2) = 0.17$
,
$ \zeta (2) = -0.27$
and
$ \dot {\zeta }(2) = 0.27$
. Solving (3.2) and (3.3), we can obtain the analytical cavity depth for a single-raindrop impact
$\varDelta =\alpha +\zeta$
.
Figure 11 compares the temporal evolution of cavity depth in cases SR, D2, D3 and D4. The cavity depths for all cases initially follow similar trends and agree well with the analytical prediction up to approximately 20 ms, beyond which notable differences emerge depending on the configuration. In case D2, cavity depth increases significantly after 21.7 ms and continues to grow until 27.6 ms, after which the growth stagnates. Following 36.8 ms, the cavity begins to retract, leading to a corresponding decrease in depth. Case D3 reaches its first peak cavity depth at 27.6 ms, followed by a slight decrease due to cavity stagnation and retraction. Upon bulge formation, there is a modest increase in depth after 35.9 ms, reaching a second peak at 41.3 ms before retracting again. In its first retraction phase between 27.6 and 35.9 ms, the cavity in case D3 retracts more slowly than that in case SR. This is attributed to the proximity and subsequent merging of its cavities. The cavity depth of the D4 case remains nearly identical to the D3 case up to 35.9 ms. Since bulge formation is weak, no second peak is observed. Note that the SR case shows an irregular increase in cavity depth after 36 ms, which results from the downward liquid jet that penetrates the cavity bottom captured by our cavity depth algorithm (Appendix A).
Comparison of cavity depth evolution for single- and two-raindrop impact cases against the analytical model.
In addition, to quantify the evolution of the morphology, we calculate the pool surface energy
$E_p$
, defined as the surface area of the largest connected fluid region in the domain multiplied by
$\sigma _0$
. As shown in figure 12, the normalised pool surface energy
$E_p(t)/E_p(0)$
remains nearly identical across all cases until approximately 9 ms. Beyond this time, case D2 exhibits a slower increase compared to the others, indicating fluid detachment from the pool associated with its rapid central film disintegration (figures 8
e–g). For case D3, a sharper decline appears at approximately 25 ms, corresponding to the onset of central film break-up (figures 9
e–g). In contrast, because case D4 shows no significant central film disintegration, its
$E_p$
follows the expected trend of cavity expansion and retraction without any notable deviation.
Comparison of the normalised pool surface energy
$E_p(t)/E_p(0)$
with time between D2, D3 and D4.
3.3. Raindrop interaction and droplet evolution
In this subsection, the droplet evolution is investigated, with a focus on how inter-raindrop interaction in cases D2, D3 and D4 shapes secondary droplet production, spatial distribution and re-merging dynamics. Our analysis of secondary droplet behaviour extends up to 15 ms, as they are no longer produced in significant numbers beyond this time. We distinguish between the production of droplets, their redistribution above the free surface, and their transport into the cavity, since these processes are affected differently by the central liquid film and cavity airflow.
3.3.1. Droplet production and spatial distribution
The temporal evolution of
$N_s^*=N_s/N_{\textit{rain}}$
, the number of secondary droplets normalised by the number of raindrops, at different sizes in cases SR, D2, D3 and D4, is shown in figure 13. This size–time distribution provides an overview of when droplets of different sizes are produced, and how long they remain in the domain. Overall, all cases demonstrate similar trends: small secondary droplets are produced early in the impact in large quantities, while larger droplets appear later in smaller quantities. This indicates that the general sequence of droplet production is similar across the single- and two-raindrop cases. The main spacing-dependent difference appears at later times, where case D2 shows a higher
$N_s^*$
for droplets with radii
$r_s$
in the range
$0.1\lt r_s\lt 0.5\,\mathrm{mm}$
during
$9\lt t\lt 15\,\mathrm{ms}$
. These secondary droplets originate from central film disintegration, which is absent in SR. As discussed previously, the central film rupture for cases D3 and D4 is not violent enough to create secondary droplets, leading to the lower
$N_s^*$
in this size range.
Temporal evolution of the normalised number of secondary droplets at different sizes for cases (a) SR, (b) D2, (c) D3 and (d) D4.
To determine whether the late-time larger-droplet population in D2 is indeed associated with central film disintegration, their spatial number density profiles along the
$x$
- and
$y$
-directions are plotted at selected times from 0.28 ms to 11.68 ms in figure 14. Here,
$n(x)$
and
$n(y)$
denote the number densities of secondary droplets in the
$x$
- and
$y$
-directions per unit distance, respectively. The
$x$
-direction profile is particularly useful because droplets produced by central film disintegration are expected to appear near the mid-plane between the two impact sites.
The behaviour of
$n(x)$
can be approximately divided into two stages, before and after the central film disintegration. At 0.28 ms, relatively uniform distributions, appearing as plateaus, are observed (figure 14
a). For cases D3 and D4, two distinct plateaus are observed, reflecting their larger initial raindrop spacing, while case D2 exhibits a single plateau due to its smaller inter-raindrop distance. A twin-peak distribution appears in all cases (figure 14
c) as droplets fall towards the two cavities. At 8.83 ms, although cases D3 and D4 maintain a two-peak distribution (figure 14
e), case D2 develops a third peak near
$x=0$
, which is attributed to central film disintegration producing additional droplets. By 11.68 ms (figure 14
g), the profiles show a prominent central peak for case D2, indicating that most of its secondary droplets in the domain at this time were produced by central film disintegration. In contrast, cases D3 and D4 still maintain a two-peak distribution, although with significantly fewer droplets overall. Therefore, the appearance of the central peak in D2 provides a spatial signature of secondary droplet production by central film disintegration, while the two-peak distributions in D3 and D4 indicate that their droplet populations remain primarily associated with the two original impact sites.
For
$n(y)$
, the initial profiles (figure 14
b) are sharply concentrated near the unperturbed free surface
$y=0$
, where the droplets are generated. By 2.56 ms (figure 14
d), the peak spreads in both directions, and droplets begin to appear at negative
$y$
values, indicating that they are falling into the cavity. The total number density continues to decline over time as droplets either exit the top of the computational domain or re-merge with the pool at the cavity bottom. This behaviour is reflected in the
$y$
-profiles, which evolve from a single peak to two peaks, with one moving in the positive direction, and the other in the negative direction, over time.
Number density profiles of secondary droplets along the
$x$
- and
$y$
-directions for D2, D3 and D4 cases at (a,b) 0.28 ms, (c,d) 2.56 ms, (e, f) 8.83 ms, and (g,h) 11.68 ms.
To quantify the spatial spread of the ejected droplets, we plot the normalised secondary droplet volume
$V_s/V_0$
(see figure 15), for cases SR, D2, D3 and D4 from
$t=2.85$
ms to
$t=11.97\ \mathrm{ms}$
. Here,
$V_s$
is the total secondary droplet volume within each radial and azimuthal bin, and
$V_0$
is the initial raindrop volume. The influence of inter-raindrop interaction becomes most apparent at later times. In case SR, crown convergence towards the impact centre causes the ligaments to merge into a single upward jet, limiting the lateral spread of airborne liquid (figures 15
e,i,m). In contrast, the interacting crowns in the two-raindrop cases suppress this radial convergence, allowing secondary droplets to be distributed over a broader region. Finally, the high-volume regions observed at later times (figures 15
n,o,p) should not be interpreted as an increase in droplet number. Instead, they are associated with larger secondary droplets carrying greater individual volume.
Spatial distributions of the normalised airborne secondary droplet volume
$V_s(r,\theta )/V_0$
. The columns from left to right correspond to cases SR, D2, D3 and D4, while the rows from top to bottom correspond to
$t = 2.85$
, 5.13, 9.12 and 11.97 ms. The maps are centred at
$(x,z)=(0,0)$
, and include only secondary droplets located above the undisturbed free surface. Each circular ring represents radial distance
$2d_0$
from the origin, and the white markers denote the initial impact centre locations.
Figure 15 Long description
Panel A: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case SR at t = 2.85 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel B: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D2 at t = 2.85 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel C: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D3 at t = 2.85 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel D: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D4 at t = 2.85 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel E: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case SR at t = 5.13 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel F: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D2 at t = 5.13 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel G: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D3 at t = 5.13 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel H: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D4 at t = 5.13 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel I: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case SR at t = 9.12 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel J: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D2 at t = 9.12 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel K: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D3 at t = 9.12 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel L: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D4 at t = 9.12 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel M: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case SR at t = 11.97 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel N: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D2 at t = 11.97 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel O: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D3 at t = 11.97 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue. Panel P: A heat map showing the spatial distribution of normalized airborne secondary droplet volume for case D4 at t = 11.97 ms. The map is centered at (x,z) = (0,0) and includes only secondary droplets located above the undisturbed free surface. Each circular ring represents a radial distance of 2d0 from the origin, and the white markers denote the initial impact center locations. The color scale ranges from 0 to 3.5 x 10^-3, with higher values indicated by red and lower values by blue.
3.3.2. Droplet re-merging dynamics
Here, we examine how inter-raindrop interaction affects the re-merging behaviour of secondary droplets, and separately plot
$N_s^*$
above and below the unperturbed free surface
$y=0$
for droplets with
$r_s\lt 0.05\,\mathrm{mm}$
and
$0.05\lt r_s\lt 0.1\,\mathrm{mm}$
in figure 16. Droplets below
$y=0$
are interpreted as droplets that have entered the cavity, and the decrease of their population is used to assess their re-merging with the pool.
For the droplets with
$r_s \lt 0.05\,\mathrm{mm}$
, their numbers remain largely unchanged in case SR until 15 ms (figure 13
a), whereas they gradually vanish in cases D2 and D3 (figures 13
b,c). For case D4, the number follows a trend similar to case SR, but with slightly fewer droplets overall (figure 13
d). As shown in figures 16(a) and 16(c), secondary droplets appear below
$y=0$
as early as 1 ms, consistent with their early production. Relative to case SR, the two-raindrop cases show fewer droplets remaining above
$y=0$
, and more below
$y=0$
, indicating that small droplets are more readily transported into the cavity in the two-raindrop cases, suggesting that the inter-raindrop interactions promote their downward transport.
For droplets with
$0.05\lt r_s\lt 0.1\,\mathrm{mm}$
,
$N_s^*$
decreases drastically near 9 ms for case SR (figure 13
a), a feature absent in the two-raindrop cases. Their values above and below
$y=0$
are plotted in figures 16(b) and 16(d), respectively. As shown in figure 16(d), droplets appear below the surface later than for
$r_s\lt 0.05\,\mathrm{mm}$
(figure 16
c) at approximately 4 ms, reflecting their later production. After approximately 8 ms, the values of
$N_s$
below the surface begin to diverge, with a rapid decline first in case SR, followed by cases D4, D3 and D2. Hence a similar number of secondary droplets enter the cavity across cases, but they re-merge with the pool at different rates. The order of this decline shows that the re-merging of these droplets becomes progressively slower as the inter-raindrop spacing decreases.
Time histories of the normalised number of droplets
$N_s^*$
(a,b) above and (c,d) below
$y=0$
. Here, the droplet radii are (a,c)
$r_s \lt 0.05$
mm and (b,d)
$0.05\lt r_s\lt 0.1$
mm.
The normalised vertical air velocity field
$u_{a,y}/U_0$
and normalised vorticity field
$\omega _z d_0/U_0$
are first examined inside the cavities to clarify how inter-raindrop interaction shapes the re-merging dynamics discussed above. In case SR (figures 17
a,c,e), crown formation is accompanied by vortex formation near the crown rim as the inflowing air detaches, which induces a low-pressure region and a converging crown that further strengthens the downward airflow. This qualitatively agrees with the mechanism discussed by Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2025). In contrast, in case D2 (figures 17
b,d, f), the central film blocks air from entering the inner side of each cavity, reducing the development of the downward airflow.
We then compute the normalised vertical drag experienced by the droplets,
$F_{D,y}/F_g$
, where
$F_g=(\rho _w-\rho _a)V_sg$
is the gravitational force, with
$V_s=(4/3)\pi r_s^3$
being the effective droplet volume. The vertical drag is given by
$F_{D,y}=-(\pi r_s^2/2)\rho _a C_D\, |\boldsymbol{v}_s-\boldsymbol{u}_a|\,(v_{s,y}-u_{a,y})$
, where
$\boldsymbol{v}_s$
is the droplet velocity,
$\boldsymbol{u}_a$
is the local air velocity, and
$C_D$
is estimated using the Schiller–Naumann correlation (Schiller & Naumann Reference Santini, Fest-Santini and Cossali1933). Positive values of
$F_{D,y}/F_g$
indicate upward drag from the surrounding air, whereas negative values indicate downward drag.
Normalised vertical velocity and vorticity fields at (a,b)
$t=2.28\,\mathrm{ms}$
, (c,d)
$t=4.56\,\mathrm{ms}$
and (e, f)
$t=7.41\,\mathrm{ms}$
for cases (a,c,e) SR and (b,d, f) D2 in the vertical plane
$z=0$
. The left-hand half of each plot shows the normalised vertical air velocity, while the right-hand half shows the normalised vorticity. The positions of secondary droplets with
$0.05\lt r_s\lt 0.1\,\mathrm{mm}$
are overlaid and coloured by their normalised vertical drag. The air–water interface is denoted by grey lines.
Figure 18 shows the time history of the median
$F_{D,y}/F_g$
for droplets of different sizes. For smaller droplets with
$r_s\lt 0.05\,\mathrm{mm}$
, the median
$F_{D,y}/F_g$
has large negative values at early times (figure 18
a), which is a consequence of their rapid upward ejection and small gravitational force. Later, the drag force becomes positive or nearly neutral in cases D2 and D3, as the droplets are captured into the cavity, whereas in cases SR and D4, the force remains mostly negative. This is consistent with our previous observations that a larger fraction of droplets remains outside the cavity in these two cases (figures 16
a,c).
For larger droplets with
$0.05\lt r_s\lt 0.1\,\mathrm{mm}$
, we examine those located below the undisturbed free surface (
$y=0$
) (figure 18
b). The median drag force increases much more abruptly in case SR than in cases D2 and D3, while case D4 exhibits a delayed drag increase. This behaviour can be explained by the drag forces acting on these larger droplets, as plotted in figure 17. In case SR, most droplets have fallen below
$y=0$
by
$7.4\,\mathrm{ms}$
, and no additional droplets can enter the cavity after its closure (figure 17
e). Consequently, the median drag force increases sharply as the droplets decelerate while approaching the cavity bottom. In case D2, however, the capture of droplets by the cavity (figure 17
f) continues for a much longer duration, resulting in the relatively steady median drag force observed in figure 18(b). We note that this overall trend is also consistent with the remerging behaviour discussed in figure 16(d).
Time histories of the median normalised vertical drag
$F_{D,y}/F_g$
for secondary droplets with (a)
$r_s\lt 0.05\,\mathrm{mm}$
over the entire domain, and (b)
$0.05\lt r_s\lt 0.1\,\mathrm{mm}$
located below the undisturbed free surface
$y=0$
.
4. Conclusion
In this study, the fate of secondary droplets produced by the impact of raindrops interacting with a deep liquid pool has been examined through multi-phase direct numerical simulations, accounting for varying surface tension, raindrop diameter and raindrop interactions. At an early stage of the impact, the secondary droplet size distribution has been found to adhere to the scaling law
$N_d(r_s)\propto Q(\sigma /\rho _w)^{-1/2}r_s^{-5/2}$
. Simulation data confirm this new scaling law, as well as the number density dependence on the surface tension and raindrop diameters.
The raindrop interaction has been examined to understand the spatial and temporal distributions of secondary droplets. It is found that the central film formed by the merging of adjacent rising crowns plays a key role in shaping the overall impact morphology. As the inter-raindrop distance increases, the film forms later and ruptures less violently, leading to weaker disintegration, reduced bulge formation and a shallower cavity. The impact of raindrop interaction on droplet distribution has been found to be size-sensitive: smaller droplets with
$r_s\lt 0.05\,\mathrm{mm}$
are increasingly ejected upwards rather than guided into the cavity as inter-raindrop distance grows, whereas those with
$0.05 \lt r_s \lt 0.1\,\mathrm{mm}$
disappear more quickly into the pool in the single-raindrop case than in the two-raindrop cases. The aerodynamic forcing analysis shows that this difference is associated with the cavity airflow, as the central film in the two-raindrop cases weakens the downward airflow and reduces the downward acceleration of droplets.
Having examined the secondary droplet dynamics from high-speed raindrops in detail, and established their behaviour during the interactions between adjacent raindrops, the present work provides a basis for several future directions in modelling oceanic precipitation. Higher maximum refinement levels and finer temporal resolution will allow droplets from the prompt splash to be fully captured and analysed. Expanding the computational domain would enable simulations to track the maximum travel distances of secondary droplets, providing additional insight into their influence on air–sea interactions. In addition, although the present results provide a statistical scaling for the secondary droplet population, distinguishing the respective roles of crown rim instability, ligament formation, capillary break-up and unsteady airflow remains an important direction for future work. The agreement of the two-raindrop cases with the proposed scaling suggests that the scaling may remain relevant for more general multi-raindrop impacts, at least locally where neighbouring drop interactions dominate, though its applicability beyond the single- and two-raindrop configurations requires rigorous verification in future studies. A natural extension of the simulations presented here would therefore involve a larger number of adjacent raindrops, similar to the laboratory rainfall experiments (e.g. Liu et al. Reference Liang, Zhang, Yu, Chen and Shen2024). We note, however, that this will depend on continued advancements in computing power.
Supplementary materials and movies
Supplementary materials and movies are available at https://doi.org/10.1017/jfm.2026.12028.
Acknowledgements
The anonymous referees are gratefully acknowledged for their comments and suggestions.
Funding
X.H. is supported by UC San Diego.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Impact geometry extraction
A.1. Automatic cavity geometry extraction
To quantify the cavity radius, the coordinates of interface cells located near the undisturbed free surface (
$|y| \lt 0.041$
mm) are extracted at each time step. An initial least squares circular fit is applied to these interface points in the
$x{-}z$
plane to estimate a preliminary centre and estimated radius
$r_e$
. The radial distance of each point from this estimated centre is denoted by
$r$
. To remove outliers, only points satisfying
$|r - r_e| \lt 2 \times \mathrm{median}(|r - r_e|)$
are retained. Next, a circular arc is fitted for each point using its five nearest neighbours in the
$x{-}z$
plane to estimate the local curvature. Points with curvature exceeding a predefined threshold are excluded from the fit. After these filtering steps, a final least squares circular fit is performed on the remaining points. The radius of this final circle is taken as the cavity radius
$r_c$
. The thresholds for radial deviation, curvature and neighbour count are empirically selected by iteratively adjusting them to achieve good agreement between the fitted circle and the cavity interface. An example plot at
$t = 24.51\ \mathrm{ms}$
is shown in figure 19(a).
The cavity depth is determined by first identifying interface cells in the domain, and computing their curvature. Only those with negative curvature and magnitude below a threshold
$0.122\ \mathrm{mm}^{-1}$
are retained. The five vertical neighbouring cells above each candidate are then examined; if all five are fully air, then the candidate is retained; otherwise, it is discarded. These two criteria prevent the algorithm from incorrectly identifying interface cells belonging to isolated air bubbles beneath the cavity. The minimum
$y$
-coordinate among the remaining candidates is then taken as the cavity depth. An example plot at
$t = 24.51\ \mathrm{ms}$
is shown in figure 19(b).
Example of (a) the circular fitting used to determine cavity radius, and (b) the identified cavity depth location relative to the interface point at 24.51 ms.
A.2. Automatic crown rim geometry extraction
The crown rim height and radius are determined using two different procedures because the crown geometry changes during the impact process. In the early stage, the individual ligaments around the rim remain distinct. The interface points around the crown rim are first used to extract the azimuthal variation of the rim profile (see an example in figure 20
a). The crown rim height is then defined as the median height of the locations of all ligament bases. Figure 20(b) shows the detected rim height relative to the air–water interface points. At the rim height, the radial distances of the interface points from the crown centre are computed in the
$x{-}z$
plane to estimate the rim radius.
In the later stage, after approximately
$11\ \mathrm{ms}$
, the rim converges radially, and the ligaments merge into a single liquid jet. The radial width of the interface is first extracted as a function of height, and the rim height is selected at the transition between the crown and the jet-like structure (figure 20
c). The rim radius is obtained directly from the radial-width profile at the rim height. Figure 20(d) shows the detected rim height, the air–water interface points, and the jet-like structure.
(a) Azimuthal profile of the crown rim with ligament peaks and bases, and (b) air–water interface points at
$t = 2.85\ \mathrm{ms}$
. (c) Radial width of the jet and rim, and (d) interface points at
$t = 17.67\ \mathrm{ms}$
. Also plotted are the detected rim heights.
Appendix B. Characteristics of crown corrugation and stability analysis
Following Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010), we test the feasibility of the Rayleigh–Plateau instability in predicting the growth of perturbations on the crown rim and the subsequent production of secondary droplets from these perturbations. In this framework, the local rim radius is written as
\begin{equation} r(\theta ,t)=r_0(t)+\epsilon (\theta ,t), \end{equation}
where
$r_0(t)$
is the mean rim radius,
$\theta$
is the azimuthal angle, and
$\epsilon (\theta ,t)$
is a small perturbation. The initial perturbation is decomposed into azimuthal Fourier modes,
\begin{equation} \epsilon (\theta ,0)=\sum _{m=-N}^{N} a_m\,{\rm e}^{{\rm i}m\theta }, \end{equation}
where
$m$
is the azimuthal mode number, and
$a_m$
is the initial amplitude of each mode. The physical wavenumber of mode
$m$
is
\begin{equation} k_m(t)=\frac {m}{R(t)}, \end{equation}
where
$R(t)$
is the crown radius. The dimensionless wavenumber entering the Rayleigh–Plateau dispersion relation is therefore
\begin{equation} x_m(t)=k_m(t)\,r_0(t)=\frac {m\, r_0(t)}{R(t)}. \end{equation}
For each mode, the viscous Rayleigh–Plateau growth rate
$\sigma _m$
is obtained from the classical stability result for a viscous cylindrical thread,
\begin{equation} \begin{aligned} &(2x^2)(x^2+y^2)\frac {I_1'(x)}{I_0(x)} \left [ 1-\frac {2xy}{x^2+y^2} \frac {I_1(x)}{I_1(y)} \frac {I_1'(y)}{I_1'(x)} \right ] -x^4+y^4 \\[3pt] &\qquad = \frac {\gamma r_0}{\rho \nu ^2} \frac {x\,I_1(x)}{I_0(x)}(1-x^2), \end{aligned} \end{equation}
where
$y^2=x^2+{\sigma _m r_0^2}/{\nu }$
,
$\gamma$
,
$\rho$
and
$\nu$
are the liquid surface tension, density and kinematic viscosity, respectively, and
$I_n\ (n=0,1)$
denotes the modified Bessel function of the first kind (Chandrasekhar Reference Chandrasekhar2013).
(a) Air–water interface of the crown in case SR at
$y=2.25$
mm above the unperturbed free surface at
$1.14$
ms. Here, the orange outline represents the extracted corrugation profile around the rim. (b) Power spectrum of the rim corrugation profile, i.e. the orange curve in (a), normalised by its maximum value. The vertical black dotted line marks the theoretical wavenumber with the maximum growth rate predicted by the Rayleigh–Plateau instability model.
Because the crown rim expands in time, the perturbation amplitude is affected by both capillary growth and rim stretching. Following Eggers & Villermaux (Reference Eggers and Villermaux2008), the amplitude of each mode evolves as
\begin{equation} \frac {{\rm d}\ln \epsilon _m}{{\rm d}t} = -\frac {s}{2} + \sigma _m, \quad s=\frac {\dot R}{R}. \end{equation}
After integration, this gives
\begin{equation} \epsilon _m(t) = a_m \left (\frac {R(0)}{R(t)}\right )^{1/2} \exp \left [ \int _0^t \sigma _m(t')\,{\rm d}t' \right ]\!. \end{equation}
The theoretical rim profile is then reconstructed as
\begin{equation} r(\theta ,t) = r_0(t) + \left (\frac {R(0)}{R(t)}\right )^{1/2} \sum _{m=-N}^{N} a_m \exp \left [{\rm i}m\theta + \int _0^t \sigma _m(t')\,{\rm d}t' \right ]\!. \end{equation}
Using (B8), we compute the theoretical corrugation profile and its wavenumber spectrum with
$R(t)$
and
$r_0(t)$
extracted from the simulations. The theoretical dominant wavenumber is then identified from the mode at the spectrum peak. The measured spectrum is obtained independently from the simulated interface. The interface points near the crown are projected onto the
$x$
–
$z$
plane to reconstruct the rim outline. The deviation of the extracted outer rim radius from its mean defines the measured corrugation profile
$h$
along the rim (figure 21
a) and the corresponding power spectrum (figure 21
b).
As shown, the peak wavenumber obtained from the simulation does not coincide with the theoretical prediction by the Rayleigh–Plateau model. Similar comparisons from
$t=1.14\ \mathrm{ms}$
to
$t=4\ \mathrm{ms}$
show persistent deviations (results not shown). This discrepancy is likely related to the difference between the present impact regime and the crown-splash regime studied by Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010). Their experiments used silicone oil and moderate impact velocities, giving relatively small
$\textit{Re}\sim 1000$
and
$\textit{We}\sim 800$
, and the crown corrugations evolved smoothly from small-amplitude perturbations. In contrast, our simulations occur at much higher Reynolds and Weber numbers, i.e.
$\textit{Re}\sim 30\,000$
and
$\textit{We}\sim 3000$
, where the rim deforms violently, and the crown break-up is highly fragmented and irregular (figure 21
a). Such high irregularity in the crown corrugation continues to be observed in all of our cases. Therefore, it is inappropriate to treat droplet production solely as the growth of rim corrugations through the Rayleigh–Plateau instability in the present impact regime.
Nrain
U0
dh
dv
σ/σ0
σ0
D
x
y
Lmax=10
Lmax=11
12
Nrain
(d0/d)2
Ep(t)/Ep(0)
x
y
Vs(r,θ)/V0
t=2.85
(x,z)=(0,0)
2d0
Ns∗
y=0
rs<0.05
0.05
t=2.28ms
t=4.56ms
t=7.41ms
z=0
0.05
FD,y/Fg
rs<0.05mm
0.05
y=0
t=2.85 ms
t=17.67 ms
y=2.25
1.14
Open access
