1. Introduction
When a raindrop splashes on the surface of a pond, it takes less than the blink of an eye for a crater to form beneath the surface, throwing a fluid crown into the air, and for it to collapse, propelling upwards a fluid jet. These are the key features of the splashing regime, which occurs within a specific range of drop radius, impact velocity, impact angle and physical properties of the fluids such as surface tension, density and viscosity (Rein Reference Rein1993). Worthington (Reference Worthington1908) was the first to report these features using pioneering high-speed photography methods. The splashing regime was then extensively investigated, regarding, in particular, the time evolution of the transient crater following the impact (e.g. Engel Reference Engel1967; Morton, Rudman & Jong-Leng Reference Morton, Rudman and Jong-Leng2000; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010) and the scaling of the maximum crater radius (e.g. Engel Reference Engel1966; Macklin & Metaxas Reference Macklin and Metaxas1976; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). The formation, evolution and fragmentation of the fluid crown (e.g. Allen Reference Allen1975; Krechetnikov & Homsy Reference Krechetnikov and Homsy2009; Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010; Agbaglah, Josserand & Zaleski Reference Agbaglah, Josserand and Zaleski2013) and of the central jet (e.g. Fedorchenko & Wang Reference Fedorchenko and Wang2004; Ray, Biswas & Sharma Reference Ray, Biswas and Sharma2015; van Rijn et al. Reference van Rijn, Westerweel, van Brummen, Antkowiak and Bonn2021) have also been examined.
Drop impact processes cover a wide variety of applications. These include engineering applications such as the water entry of projectiles (Clanet, Hersen & Bocquet Reference Clanet, Hersen and Bocquet2004) or spray painting (Hines Reference Hines1966). They also include Earth sciences applications such as the production of oily marine aerosol by raindrops (Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015), spray generation from raindrop impacts on seawater and soil (Zhou et al. Reference Zhou, Wang, Lu, Chen, Wang, Chen, Yang and Wang2020) and planetary impact craters (Melosh Reference Melosh1989; Landeau et al. Reference Landeau, Deguen, Phillips, Neufeld, Lherm and Dalziel2021; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2021, Reference Lherm, Deguen, Alboussière and Landeau2022). Planetary impacts occur on terrestrial planets from the early stages of accretion to modern meteorite impacts. During planetary formation, thermal and chemical partitioning between the core and the mantle is influenced by the physical mechanisms of segregation between the metal of the impactors’ cores and the silicates of the growing planet (Stevenson Reference Stevenson1990; Rubie, Nimmo & Melosh Reference Rubie, Nimmo and Melosh2015; Lherm & Deguen Reference Lherm and Deguen2018), with major implications for the chemical, thermal and magnetic evolution of the planet (Fischer et al. Reference Fischer, Nakajima, Campbell, Frost, Harries, Langenhorst, Miyajima, Pollok and Rubie2015; Badro et al. Reference Badro, Aubert, Hirose, Nomura, Blanchard, Borensztajn and Siebert2018; Olson, Sharp & Garai Reference Olson, Sharp and Garai2022). In particular, the cratering process is responsible for the initial dispersion and mixing of the impactors’ cores (Landeau et al. Reference Landeau, Deguen, Phillips, Neufeld, Lherm and Dalziel2021; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). In planetary science, impact craters are also a tool for sampling the shallow interior of planets and satellites by combining observations of planetary surfaces with excavation and ejecta deposition models (Maxwell Reference Maxwell1977; Barnhart & Nimmo Reference Barnhart and Nimmo2011; Kurosawa & Takada Reference Kurosawa and Takada2019). Therefore, understanding the implications of these planetary impacts requires the modelling of the velocity field produced during the formation of the craters.
In the splashing regime, the fate of the crater, the fluid crown and the central jet is directly related to the velocity field produced around the crater boundary. The dynamics of the crater is indeed closely related to the velocity field in the ambient fluid, in particular regarding the evolution of the shape of the cavity. The formation of the fluid crown is also related to the ambient velocity field through the mass flux distribution across the initial water surface. Finally, the production of the central jet is associated with a convergent velocity field, resulting from the collapse of the crater due in part to buoyancy forces.
The velocity field associated with crater evolution in the splashing regime has been investigated both experimentally and numerically in previous studies. Engel (Reference Engel1962) was the first to examine the velocity field around a crater by seeding the flow with particles in order to visualise the flow streamlines. These observations allowed the determination of the velocity field configuration associated with the crater expansion and its subsequent collapse. More recently, the velocity field was investigated using modern particle image velocimetry (PIV) methods. These velocity field measurements have been used to investigate the origin of vortex rings beneath a crater (Liow & Cole Reference Liow and Cole2009), the formation of the central jet (van Rijn et al. Reference van Rijn, Westerweel, van Brummen, Antkowiak and Bonn2021) or solutocapillary flows following the impact of drops on salted water (Musunuri et al. Reference Musunuri, Benouaguef, Amah, Blackmore, Fischer and Singh2017). Numerical simulations have also focused on the crater velocity field, regarding in particular the entrapment of air bubbles when the crater collapses and the formation of the central jet (Morton et al. Reference Morton, Rudman and Jong-Leng2000; Ray et al. Reference Ray, Biswas and Sharma2015).
Most of the models involving a prediction of the crater velocity field assume either an arbitrary velocity field (Maxwell Reference Maxwell1977) or an arbitrary velocity potential associated with an imposed crater geometry, such as a hemispherical crater (Engel Reference Engel1967; Leng Reference Leng2001) or a spherical crater able to translate vertically (Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010). Since these models have only been compared with experimental measurements of crater size and/or shape, a comparison with experimental measurements of the velocity field is thus required to assess their accuracy. In any event, a new model is required to consistently model the geometry of the cavity without the use of an arbitrarily imposed velocity field or potential.
In this paper, we examine simultaneously the dynamics of the cavity and of the velocity field produced in drop impact experiments. In § 2, we present the experimental set-up, methods and diagnostics, as well as the set of dimensionless numbers used in this study. In § 3, we describe the experimental results obtained for the crater shape and the velocity field. In § 4, we compare the existing velocity field models with our experimental measurements. In § 5, we finally derive a Legendre polynomial model based on an unsteady Bernoulli equation coupled with a kinematic boundary condition.
2. Experiments
2.1. Experimental set-up
In these experiments, we release a liquid drop in the air above a deep liquid pool of the same liquid (figure 1). We vary the impact velocity 
$U_i$ by changing the release height of the drop while keeping the drop radius 
$R_i$ fixed. We also keep constant the density 
$\rho$, the viscosity 
$\mu$ and the surface tension 
$\sigma$ of the fluids.
Figure 1. Schematic view of the drop impact experimental set-up.
The liquid pool is contained in a 
$16\,{\rm cm} \times 16\,{\rm cm} \times 30$ cm glass tank. The pool level is set at the top of the tank to minimise the thickness of the meniscus on the sides of the tank. This allows the imaging of a field of view unperturbed by the free-surface meniscus effect. We generate the drops using a needle supplied with fluid by a syringe driver. When the weight of the drop exceeds the surface tension forces, the drop comes off. We use a nylon plastic needle with an inner diameter of 4.7 mm, generating drops with a radius 
$R_i=2.7\,\mathrm {mm}$. We measure the drop size based on a calibration using mass measurements of dozens of drops and assuming each drop is spherical. We validate this method using high-speed pictures of the drop prior to impact where we can directly measure the drop radius. We obtain a relative difference of 1.4 % between mass measurements and direct measurements. Impact velocities are in the range 
$U_i=1\unicode{x2013}5\,\mathrm {m}\,\mathrm {s}^{-1}$. We calculate the impact velocity for each experiment using a calibrated free-fall model for the drop, including a quadratic drag. We validate this method using high-speed pictures of the drop prior to impact where we can directly measure the drop velocity. We obtain a relative difference of 0.6 % between the velocity model and direct measurements. We use water both in the drop and in the pool, in a temperature-controlled environment. The density is 
$\rho =998\pm 1\,\mathrm {kg}\,\mathrm {m}^{-3}$, measured using an Anton Paar DMA 35 Basic densitometer. The viscosity is 
$\mu =1\pm 0.01\,\mathrm {mPa}\,\mathrm {s}$ (Haynes Reference Haynes2016). The surface tension at the air–water interface is 
${\sigma =72.8\pm 0.4\,\mathrm{mJ}\,\mathrm{m}^{-2}}$ (Haynes Reference Haynes2016).
In our experiments, we position the camera at the same height as the water surface. We record images at 1400 Hz with a 
$2560\,{\rm pixel}\times 1600$ pixel resolution (
$21\,\mathrm {\mu }\mathrm {m}\,\mathrm {px}^{-1}$) and a 12 bits dynamic range, using a high-speed Phantom VEO 640L camera and a Tokina AT-X M100 PRO D Macro lens.
2.2. Dimensionless numbers
In these experiments, the impact dynamics depends on 
$U_i$, 
$R_i$, 
$\rho$, 
$\mu$, 
$\sigma$ and the acceleration of gravity 
$g$. Since these six parameters contain three fundamental units, the Vaschy–Buckingham theorem dictates that the impact dynamics depends on a set of three independent dimensionless numbers. We choose the following set:
\begin{equation} Fr=\frac{U_i^2}{g R_i}, \quad We=\frac{\rho U_i^2R_i}{\sigma}, \quad Re=\frac{\rho U_i R_i}{\mu}. \end{equation}
The Froude number 
$Fr$ is a measure of the relative importance of impactor inertia and gravity forces. It can also be interpreted as the ratio of the kinetic energy 
$\rho R_i^3 U_i^2$ of the impactor to its gravitational potential energy 
$\rho g R_i^4$ just before impact. The Weber number 
$We$ compares the impactor inertia and interfacial tension at the air–liquid interface. The Reynolds number 
$Re$ is the ratio between inertial and viscous forces. In the following, time, lengths and velocities are made dimensionless using the drop radius and the impact velocity, i.e. using respectively 
$R_i/U_i$, 
$R_i$ and 
$U_i$. These dimensionless quantities are denoted with a tilde. For example, we use a dimensionless time 
$\tilde {t}=t/(R_i/U_i)$.
We focus on four cases with Froude numbers, Weber numbers and Reynolds numbers respectively in the range 
$Fr \simeq 100\unicode{x2013}1000$, 
$We \simeq 100\unicode{x2013}1000$ and 
$Re \simeq 4400\unicode{x2013}13\,600$ (table 1). For each case, we conducted three acquisitions, with similar experimental results regarding both the crater shape (e.g. figure 4) and the velocity field (e.g. figure 11). This validates the repeatability of the experiments.
Table 1. Dimensionless numbers used in the experiments (see (2.1a–c) for details).
2.3. Particle image velocimetry
The velocity field is obtained using PIV. We seed the tank with polyamide particles (figure 1), the concentration, diameter and density of which being respectively 
$C_p=0.26\,\mathrm {g}\,\mathrm {L}^{-1}$, 
$d_p=20\,\mathrm {\mu }\mathrm {m}$ and 
$\rho _p=1030\,\mathrm {kg}\,\mathrm {m}^{-3}$. We illuminate these particles in suspension with a 
$1\,\mathrm {mm}$ thick laser sheet (
$532\,\mathrm {nm}$), produced using a continuous 
$10\,\mathrm {W}$ Nd:YAG laser, together with a diverging cylindrical lens and a telescope. The laser sheet is verticalised using a 
$45^\circ$ inclined mirror located below the tank. The laser wavelength is isolated using a band-pass filter (
$532 \pm 10\,\mathrm {nm}$).
In order to calculate the velocity field, the camera records two images of the field of view separated by a short time (
$\Delta t=200\,\mathrm {\mu }\mathrm {s}$). These two images are divided into interrogation windows in which a cross-correlation operation allows obtaining the average particle displacement. This involves a five-stage multi-pass processing with interrogation windows decreasing in size. The final interrogation window size is a 64 px square with an overlap of 75 %. In each window, a velocity vector is then calculated, which allows the construction of the velocity field over the whole field of view. Finally, the velocity field is spatially calibrated using a sight.
2.4. Experimental diagnostic
2.4.1. Crater shape
The crater shape is directly obtained from the raw images used in the PIV procedure (figure 2). The crater corresponds to a particle-free area, together with a high-light-intensity area, explained by reflections at the air–water interface, in particular at the bottom of the crater. The crater boundary is defined using these image properties, which allow the delineation of the cavity using background removal, an intensity threshold method and image binarisation.
Figure 2. Velocity field resulting from the PIV procedure, superposed on a corresponding experimental raw image. The solid lines correspond to the shape of the crater obtained from the shifted Legendre polynomial decomposition (2.2), for degrees 
$k=0$ (blue), 
$k=1$ (orange) and 
$k=2$ (green). The definitions of the Cartesian 
$(x,y,z)$ (black) and of the spherical 
$(r,\theta,\varphi )$ (red) coordinate systems are also represented.
We fit the crater boundary position 
$R$ (figure 2), which depends on the polar angle 
$\theta$ and time 
$t$, using a set of shifted Legendre polynomials 
$\bar {P}_k$ up to degree 
$k_{max}=2$:
\begin{equation} R(\theta,t)=\sum_{k=0}^{k_{max}} R_k(t)\bar{P}_k(\cos\theta), \end{equation}
where 
$R_k(t)$ are coefficients fitted with a least-squares method. The shifted Legendre polynomials are an affine transformation of the standard Legendre polynomials 
$\bar {P}_k(x)=P_k(2x-1)$, and are orthogonal on 
$[0,1]$, i.e. on a half-space. The coefficients 
$R_k(t)$ correspond to increasingly small-scale deviations from a hemispherical shape. Coefficient 
$R_0(t)$ corresponds to the mean crater radius (figure 2, blue line). Coefficient 
$R_1(t)$ corresponds to a deformation of the crater, linear in 
$\cos \theta$, with respect to an hemisphere (figure 2, orange line). When 
$R_1(t)>0$, the crater is stretched vertically, resulting in a prolate cavity. When 
$R_1(t)<0$, the crater is stretched horizontally, resulting in an oblate cavity. Finally, 
$R_2(t)$ corresponds to a deformation of the crater, quadratic in 
$\cos \theta$, with respect to a hemisphere (figure 2, green line).
In order to validate the crater shape determination procedure, we compare the coefficients 
$R_k(t)$ obtained from the raw images used in the PIV procedure (e.g. figure 2) with the coefficients obtained from an experiment under the same conditions, but illuminated from behind (e.g. figure 13). This backlight experiment (see Lherm et al. (Reference Lherm, Deguen, Alboussière and Landeau2022) for experimental details) allows reliable determination of the shape of the crater. Figure 3 shows that the coefficients are very similar between the two methods, which validates the crater shape determination procedure from PIV raw images.
Figure 3. Coefficients 
$R_0(t)/R_i$ (a), 
$R_1(t)/R_0$ (b) and 
$R_2(t)/R_0$ (c) as a function of time 
$\tilde {t}$. The circles and the solid lines correspond to the crater shape obtained respectively from the PIV procedure (case B, 
$Fr=444$) and a similar backlight experiment (
$Fr=442$).
2.4.2. Velocity field
We aim to compare the experimental velocity field obtained using the PIV procedure with velocity models. For that purpose, the velocity field 
$\boldsymbol {u}=(u_r,u_\theta,u_\varphi )$ is expressed in a spherical coordinate system (
$r,\theta,\varphi$) defined such that 
$u_r$ and 
$u_\theta$ are in the plane of the laser sheet (figure 2, red coordinates). The origin of this coordinate system is the contact point between the impacting drop and the target liquid (figure 2, point 
$O$).
We decompose the components of the velocity field on a shifted Legendre polynomial basis:
\begin{align} u_r(r, \theta, t)&=\sum_{l=0}^{+\infty} u_{r,l}(r, t) \bar{P}_l(\cos\theta), \end{align}
\begin{align} u_\theta(r, \theta, t)&=\sum_{l=0}^{+\infty} u_{\theta,l}(r, t) \bar{P}_l(\cos\theta), \end{align}
where 
$u_{r,l}(r,t)$ and 
$u_{\theta,l}(r,t)$ are respectively the decomposition coefficients of 
$u_r$ and 
$u_\theta$. The shifted Legendre polynomials 
$\bar {P}_l(\cos \theta )$ being orthogonal on half-hemispheres (
$\theta \in [0, {\rm \pi}/2]$), we obtain the 
$u_{r,l}(r,t)$ and 
$u_{\theta,l}(r,t)$ coefficients using a least-squares inversion of the experimental velocity components over the separate half-hemispheres 
$\theta \geq 0$ and 
$\theta < 0$, before averaging the results from the left and right half-hemispheres. Since the flow is close to axisymmetric (e.g. figure 2), the coefficients obtained by the inversion over each half-hemisphere are very close to each other. Assuming an axisymmetric flow, note that 
$u_{r,0}(r,t)$ is the average of 
$u_r$ over the full hemisphere.
3. Experimental results
3.1. Crater shape
Figure 4 shows the fitted coefficients of the shifted Legendre decomposition of the crater boundary (2.2) as a function of time, for all experimental cases. We normalise the fitted coefficients 
$R_1(t)$ and 
$R_2(t)$ by 
$R_0(t)$, i.e. the mean crater radius. Using this normalisation, we quantify the deviation of the crater geometry from a hemisphere. We also normalise time by the opening time scale of the crater (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022):
\begin{equation} \tilde{t}_{max}=\frac{1}{2}\left(\frac{8}{3}\right)^{1/8}\mathrm{B}\left(\frac{1}{2},\frac{5}{8}\right) \varPhi^{1/8} \xi^{1/2} Fr^{5/8}, \end{equation}
where 
$\varPhi$ and 
$\xi$ are respectively energy partitioning and kinetic energy correction coefficients, and 
$\mathrm {B}$ is the beta function. This scaling is obtained by using an energy conservation equation where the sum of the potential energy of the crater and of the kinetic energy of the crater, corrected by 
$\xi$, is equal at any instant of time to the kinetic energy of the impacting drop, corrected by 
$\varPhi$. Assuming that the kinetic energy of the crater vanishes when the cavity reaches its maximum size (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022), the maximum crater radius scales as
\begin{equation} \tilde{R}_{max}=\left(\tfrac{8}{3}\right)^{1/4}\varPhi^{1/4} Fr^{1/4}. \end{equation}
Using this 
$Fr^{1/4}$ scaling law, the energy conservation equation is integrated between 
$\tilde {t}=0$ and 
$\tilde {t}=\tilde {t}_{max}$ to obtain the opening time scale of the crater given by (3.1). More details can be found in Lherm et al. (Reference Lherm, Deguen, Alboussière and Landeau2022). With our experimental range of Froude number, we use 
$\varPhi =Fr^{-0.156}$ and 
$\xi =0.34$ (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). This normalisation allows us to collapse our experiments on the same time scale.
Figure 4. Coefficients 
$R_0(t)/R_i$ (a), 
$R_1(t)/R_0(t)$ (b) and 
$R_2(t)/R_0(t)$ (c) as a function of time normalised by the opening time scale of the crater 
$\tilde {t}/\tilde {t}_{max}$, in the four cases. Inset: maximum mean crater radius 
$R_{0_{max}}/R_i$ as a function of the Froude number 
$Fr$. For each case, the different types of markers correspond to different experiments.
In figure 4, the crater shape evolution of case A is markedly different from that of cases B, C and D. Thus, we describe this case separately. We first deal with the high-
$We$ experiments (cases B, C and D), where surface tension effects are negligible in comparison with the impactor inertia (e.g. Pumphrey & Elmore Reference Pumphrey and Elmore1990; Morton et al. Reference Morton, Rudman and Jong-Leng2000; Leng Reference Leng2001; Ray et al. Reference Ray, Biswas and Sharma2015). The crater size increases with the Froude number (figure 4, inset) in a way that is compatible with a 
$Fr^{1/4}$ scaling law for the maximum mean crater radius 
$R_{0_{max}}$ (Engel Reference Engel1966; Leng Reference Leng2001; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). Furthermore, the evolution of the crater shape relative to the mean crater size is independent of the Froude number, with similar evolution of 
$R_1(t)/R_0(t)$ and 
$R_2(t)/R_0(t)$ (figure 4b,c).
At early times of the crater opening stage (
$\tilde {t}/\tilde {t}_{max} \lesssim 0.25$), the mean radius of the crater 
$R_0(t)$ increases (figure 4a) as the cavity opens. The crater has a flat-bottomed oblate shape (e.g. figure 5i) as a result of the spread of the drop on the surface of the pool, with negative 
$R_1(t)$ (figure 4b). The flat-bottomed oblate cavity gradually becomes hemispherical as a result of the overpressure produced at the contact point between the impacting drop and the surface (e.g. figure 5ii). The magnitude of 
$R_1(t)/R_0(t)$ indeed decreases with time during this stage (figure 4b). The crater is also deformed at higher degrees with mostly negative 
$R_2(t)/R_0(t)$ (figure 4c). This corresponds to second-order deviations from the hemispherical shape, with a flattened crater boundary close to the surface (e.g. figure 5i).
Figure 5. Time evolution of the velocity 
$|\boldsymbol {u}|=(u_x^2+u_z^2)^{1/2}$ for case B. The vector field corresponds to the experimental velocity field, normalised by its maximum value in each snapshot. The solid green line corresponds to the crater boundary determined using the Legendre polynomial decomposition (2.2).
At intermediate times of the crater opening stage (
$0.25 \lesssim \tilde {t}/\tilde {t}_{max} \lesssim 0.5$), the crater continues to open (figure 4a). The cavity is still stretched vertically, which leads to increasingly positive 
$R_1(t)/R_0(t)$ (figure 4b), i.e. a prolate cavity (figure 5iii). The crater reaches a maximum prolate deformation when 
$\tilde {t}/\tilde {t}_{max} \simeq 0.5$, with 
$R_1(t)/R_0(t) \simeq 0.08$ (figure 4b). The crater is also deformed at higher degrees, with positive 
$R_2(t)/R_0(t)$ (figure 4c). This corresponds to a vertical crater boundary close to the surface (figure 5iii).
At late times of the crater opening phase (
$0.5 \lesssim \tilde {t}/\tilde {t}_{max} \lesssim 1$), the mean crater radius still increases (figure 4a) but the crater starts to flatten with decreasing 
$R_1(t)/R_0(t)$ (figure 4b). As the opening velocity of the crater decreases, buoyancy forces become significant, resulting in the horizontal stretching of the cavity. The crater flattens to give an approximately hemispherical crater at 
$\tilde {t}/\tilde {t}_{max} \simeq 1$ (figure 5v).
After the crater has reached its maximum size (
$\tilde {t}/\tilde {t}_{max} \gtrsim 1$), the mean crater radius starts to decrease (figure 4a). The value of 
$R_1(t)/R_0(t)$ decreases at a rate higher than in the opening stage of the crater (figure 4b). Horizontal stretching of the crater is accelerated, as expected since buoyancy forces are now prevailing. This leads to the formation of an increasingly oblate cavity (figure 5vi–vii). When 
$\tilde {t}/\tilde {t}_{max} \gtrsim 1.5$, higher degrees eventually deviate from zero with positive 
$R_2(t)/R_0(t)$ (figure 4c). In addition to the negative value of 
$R_1(t)/R_0(t)$, this corresponds to the formation of the central jet (figure 5viii–ix).
We now deal with the moderate-
$We$ experiment (case A), where surface tension effects are significant in comparison with the impactor inertia (e.g. Pumphrey & Elmore Reference Pumphrey and Elmore1990; Morton et al. Reference Morton, Rudman and Jong-Leng2000; Leng Reference Leng2001; Ray et al. Reference Ray, Biswas and Sharma2015). In this case, a downward-propagating capillary wave develops at the cavity interface and drives the crater deformation, often leading to the entrapment of a bubble due to the pinching of the cavity (e.g. Oguz & Prosperetti Reference Oguz and Prosperetti1990; Pumphrey & Elmore Reference Pumphrey and Elmore1990; Prosperetti & Oguz Reference Prosperetti and Oguz1993; Elmore, Chahine & Oguz Reference Elmore, Chahine and Oguz2001). This mechanism is typically expected at moderate 
$We$, i.e. 
$We \simeq 30$–
$140$ (figure 6 in Pumphrey & Elmore Reference Pumphrey and Elmore1990). During crater opening, this explains why the maximum prolate deformation occurs later than in the other cases, at 
$\tilde {t}/\tilde {t}_{max} \simeq 0.8$, and why the prolate deformation is larger, with 
$R_1(t)/R_0(t) \simeq 0.2$ (figure 4b). During crater closing, the evolution of 
$R_1(t)/R_0(t)$ and 
$R_2(t)/R_0(t)$ is markedly different from that of the other cases due to the convergence of the capillary wave at the bottom of the crater.
3.2. Velocity field
3.2.1. Velocity maps
Figure 5 shows the evolution of the norm of the velocity 
$|\boldsymbol {u}|=(u_x^2+u_z^2)^{1/2}$ as a function of time, for case B. During the opening stage of the crater, the velocity around the crater gradually decreases due to the deceleration of the crater boundary (figure 5i–iv). The maximum velocity is 
${\simeq }1.1\,\mathrm {m}\,\mathrm {s}^{-1}$ at time 1.3 ms after contact (figure 5i), which corresponds to 32 % of the impact velocity. When 
$\tilde {t}/\tilde {t}_{max} \gtrsim 0.1$, the norm of the velocity decreases radially around the crater (figure 5ii–iv), whereas, when 
$\tilde {t}/\tilde {t}_{max} \lesssim 0.1$, the velocity decreases at a higher rate on the side of the crater. This may be explained by the initial oblate shape of the crater, related to the spread of the drop on the water surface upon impact, which leads to a higher velocity beneath the crater as it becomes gradually hemispherical. The velocity field is composed of a dominant radial component and of a polar component responsible for an upward flow across the initial water surface (figure 5i–iv). The polar component is thus responsible for the formation of the liquid crown above the water surface (e.g. Rein Reference Rein1993; Fedorchenko & Wang Reference Fedorchenko and Wang2004; Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010).
When the crater reaches its maximum size (figure 5v), the cavity is nearly hemispherical and the velocity field seems to vanish simultaneously in the entire flow, consistent with the observations of Engel (Reference Engel1966), which were subsequently used in several velocity models (e.g. Engel Reference Engel1967; Prosperetti & Oguz Reference Prosperetti and Oguz1993). However, this first-order assumption on the simultaneous vanishing velocity field does not hold when the flow is examined in detail. Beneath the cavity, the velocity gradually decreases and eventually vanishes just before the crater reaches its maximum size. The velocity is directed downwards due to the expansion of the crater. The velocity then increases again but is directed upwards due to the collapse of the crater. On the side of the cavity, close to the surface, the velocity does not vanish when the crater reaches its maximum size. The collapse of the crater takes over its initial expansion, which allows the retention of outward velocities on the side of the crater.
When the crater collapses (figure 5vi–ix), a convergent flow forms towards the centre of the cavity. This leads to the formation of the central jet.
Figure 6 shows the evolution of the vorticity 
$\omega _y=\partial u_x/\partial z-\partial u_z/\partial x$ as a function of time, for case B. The vorticity produced by the impact around the crater is confined close to the air–water boundary, in particular when the crater is strongly deformed, at the beginning of the crater opening (figure 6i–ii) and when it collapses (figure 6vi–ix). This suggests that the flow is mostly irrotational, which supports the potential flow assumption used in previous models (§ 4). Furthermore, some of the vorticity observed around the crater boundary may be an artefact related to spurious velocity measurements produced by cross-correlations on reflections at the air–water interface, and not on PIV particles. This assumption is supported by the estimated diffusion length of the vorticity (0.3 mm in 100 ms) which is significantly smaller than the typical size of the vorticity band.
Figure 6. Time evolution of the vorticity 
$\omega _y=\partial u_x/\partial z-\partial u_z/\partial x$ for case B. The vector field corresponds to the experimental velocity field, normalised by its maximum value in each snapshot. The solid green line corresponds to the crater boundary determined using the Legendre polynomial decomposition (2.2).
3.2.2. Velocity coefficients
Figure 7 shows the coefficients 
$u_{r,l}(r,t)$ and 
$u_{\theta,l}(r,t)$ (equations (2.3)–(2.4)) as a function of the radial coordinate at a given time 
$\tilde {t}=15.4$ (
$\tilde {t}/\tilde {t}_{max}=0.43$) during the crater opening stage of case B. During this stage, the velocity field is dominated by the degrees 
$l=0$ and 
$l=1$, the higher degrees 
$l\geq 2$ being much smaller. When 
$r \lesssim 1.2R_0$, we observe a decrease in the slope of the coefficients. This may be related to the deviation of the crater from a hemisphere. The coefficients indeed sample points located at varying distances from the actual crater boundary, including artefacts located into the crater, which may influence the radial dependency of these coefficients close to the crater boundary. In figure 7, we identify this misleading trend by using dashed lines when the radius is smaller than 
$\max \{R(\theta )\}$ (
$r \leq 1.11 R_0$ in figure 7).
Figure 7. Coefficients 
$u_{r,l}(r,t)$ and 
$u_{\theta,l}(r,t)$ normalised by the mean crater velocity 
$\dot {R}_0(t)$, as a function of the radial coordinate 
$r$, normalised by the mean crater radius 
$R_0(t)$, up to degree 
$l=5$. Dashed lines correspond to regions where the coefficients sample velocity artefacts located in the crater. The coefficients are calculated at 
$\tilde {t}=15.4$ (
$\tilde {t}/\tilde {t}_{max}=0.43$) for case B.
Figures 8 and 9 compare the time evolution of the coefficients 
$u_{r,l}(r,t)$ and 
$u_{\theta,l}(r,t)$ (equations (2.3)–(2.4)) between the cases, for 
$l \leq 2$. Except for the different normalisation, figure 7 is thus similar to a radial slice of these coefficient maps, for case B, at 
$\tilde {t}/\tilde {t}_{max}=0.43$. As for the crater shape, the moderate-
$We$ case A is different from the high-
$We$ cases B, C and D, for both the radial (figure 8) and the polar (figure 9) components of the velocity field. We thus deal with this case separately.
Figure 8. Time evolution of the coefficient 
$u_{r,l}(r,t)$ (
$l\in \{0,1,2\}$) normalised by the impact velocity 
$U_i$, as a function of the radial coordinate 
$r$ normalised by the drop radius 
$R_i$, for case B. Time is normalised by the opening time scale of the crater 
$\tilde {t}_{max}$.
Figure 9. Time evolution of the coefficient 
$u_{\theta,l}(r,t)$ (
$l\in \{0,1,2\}$) normalised by the impact velocity 
$U_i$, as a function of the radial coordinate 
$r$ normalised by the drop radius 
$R_i$, for case B. Time is normalised by the opening time scale of the crater 
$\tilde {t}_{max}$.
We first deal with the high-
$We$ experiments (cases B, C and D), where both components of the velocity field are similar among cases, regardless of the degree in question. The velocity field is mostly dominated by the degrees 
$l=0$ and 
$l=1$, during both the opening and the closing stage of the crater, in agreement with figure 7.
During the crater opening stage (
$\tilde {t}/\tilde {t}_{max} \lesssim 1$), the dominant degrees of the radial component 
$u_{r,0}(r,t)$ and 
$u_{r,1}(r,t)$ are positive (figure 8). This corresponds to the strong radial velocity field related to the expansion of the cavity. The dominant degrees of the polar component 
$u_{\theta,0}(r,t)$ and 
$u_{\theta,1}(r,t)$ are concomitantly positive and negative, respectively (figure 9), with a lower magnitude. This corresponds to a polar perturbation of the dominant radial velocity field, related to the mass flux across the surface 
$z=0$ which produces the fluid crown. The positive coefficient 
$u_{\theta,0}(r,t)$ indeed corresponds to a flow towards the surface, while the negative coefficient 
$u_{\theta,1}(r,t)$ corresponds to a degree 
$l=1$ perturbation, linear in 
$\cos \theta$. The degree 
$l=2$ also contributes to the velocity field of both components, in particular when the crater is strongly deformed due to the spread of the drop at the surface of the pool, at the beginning of the opening stage (
$\tilde {t}/\tilde {t}_{max} \lesssim 0.25$).
When the crater reaches its maximum size (
$\tilde {t}/\tilde {t}_{max} \simeq 1$), the dominant degrees of both components change signs as the crater starts to collapse. In detail, 
$u_{r,0}(r,t)$ vanishes later (
$\tilde {t}/\tilde {t}_{max} \simeq 1$) than 
$u_{r,1}(r,t)$ (
$\tilde {t}/\tilde {t}_{max} \simeq 0.6$) (figure 8) and 
$u_{\theta,0}(r,t)$ (
$\tilde {t}/\tilde {t}_{max} \simeq 0.8$) (figure 9). This is in agreement with the observations of figure 5 at 
$\tilde {t}/\tilde {t}_{max} \simeq 1$, where the velocity vanishes beneath the crater but not on the sides.
During the crater closing stage (
$\tilde {t}/\tilde {t}_{max} \gtrsim 1$), 
$u_{r,0}(r,t)$ and 
$u_{r,1}(r,t)$ are both negative (figure 8) and 
$u_{\theta,0}(r,t)$ is negative (figure 9). This corresponds to the development of the convergent flow related to the collapse of the crater and the formation of the central jet. As at the beginning of the opening stage, the degree 
$l=2$ of both components contributes significantly to the velocity field at the end of the closing stage (
$\tilde {t}/\tilde {t}_{max} \gtrsim 1.5$), in relation to the strongly deformed crater boundary.
We now deal with the moderate-
$We$ experiment (case A). Although the degree 
$l=0$ of both components is similar to that of the high-
$We$ cases, the degrees 
$l=1$ and 
$l=2$ of case A are significantly larger than their counterparts of cases B, C and D. Furthermore, the time at which 
$u_{r,1}(r,t)$ and 
$u_{\theta,0}(r,t)$ vanish is significantly modified. This may also be a consequence of significant surface tension effects in this moderate-
$We$ experiment, related to vigorous deformations of the crater boundary by the propagation of a capillary wave towards the bottom of the crater.
4. Comparison with existing velocity models
In this section, we review the velocity models proposed by Engel (Reference Engel1967), Maxwell (Reference Maxwell1977), Leng (Reference Leng2001) and Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010), and compare their predictions with our observations. Since most of these models have been designed to understand the crater opening stage, we compare these models with our experimental velocity measurements by focusing on a typical snapshot of this initial stage. For that purpose, figure 10 shows the dominant coefficients 
$u_{r,0}(r,t)$, 
$u_{r,1}(r,t)$, 
$u_{\theta,0}(r,t)$ and 
$u_{\theta,1}(r,t)$ of case B as a function of the radial coordinate at 
$\tilde {t}/\tilde {t}_{max}=0.24$, as well as the predictions of the models.
Figure 10. Coefficients 
$u_{r,0}(r,t)$ (a), 
$u_{r,1}(r,t)$ (b), 
$u_{\theta,0}(r,t)$ (c) and 
$u_{\theta,1}(r,t)$ (d) normalised by the mean crater velocity 
$\dot {R}_0(t)$, as a function of the radial coordinate 
$r$, normalised by the mean crater radius 
$R_0(t)$. The circles correspond to case B at 
$\tilde {t}/\tilde {t}_{max}=0.24$. The dash-dotted lines correspond to the models of Engel (Reference Engel1967), Maxwell (Reference Maxwell1977), Leng (Reference Leng2001) and Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010). The solid line corresponds to the solution of the predictive model, using the simplified equation system (5.13) with the reference set of initial conditions (5.14).
4.1. The model of Engel (1967)
The model of Engel (Reference Engel1967) assumes an energy balance where the potential energy of the crater, the potential energy of a cylindrical wave developing above the surface, the surface tension energy of the produced interface, the kinetic energy of the flow around the crater, the kinetic energy of the cylindrical wave and viscous dissipation are equal at any time to half of the kinetic energy of the impacting drop. Among the assumptions of such a model, Engel (Reference Engel1967) assumes a hemispherical crater with a radius 
$R_0(t)$ and a potential flow with a velocity potential 
$\phi$ satisfying the boundary conditions on the velocity 
$|\boldsymbol {u}|(r=+\infty )=0$ and 
$|\boldsymbol {u}|[r=R_0(t)]=\dot {R}_0(t)$. The velocity potential used in the model is
\begin{equation} \phi={-}\frac{\dot{R}_0R_0^2\cos\theta}{r}. \end{equation}
The radial component 
$u_r$ and the polar component 
$u_\theta$ of the velocity field, obtained by deriving the velocity potential, are written as
\begin{equation} \left. \begin{aligned} & u_r=\frac{\dot{R}_0R_0^2\cos\theta}{r^2},\\ & u_\theta=\frac{\dot{R}_0R_0^2\sin\theta}{r^2}. \end{aligned} \right\} \end{equation} This model allows one to capture the evolution of the mean crater radius (e.g. Engel Reference Engel1967, figure 3). The velocity field has 
$l=0$ (figure 10a) and 
$l=1$ (figure 10b) radial components and 
$l=0$ (figure 10c) and 
$l=1$ (figure 10d) polar components. This allows one to obtain a velocity field qualitatively similar to that of the experiments, including in particular a degree 
$l=1$ of the radial component, and a polar component. However, the slopes of the velocity components are smaller than the experimental slopes, in particular the 
$1/r^2$ slope of 
$u_{r,0}(r,t)$. The main limitations of the model of Engel (Reference Engel1967) are the fixed hemispherical geometry of the crater and the arbitrary velocity potential defined to fit experimental observations of the velocity field. More importantly, this velocity potential (4.1) corresponds to the flow around an expanding cylinder (with 
$r$ being the distance from the cylinder axis and 
$\theta$ the angular position around this axis) rather than around an expanding sphere, as incorrectly assumed in Engel (Reference Engel1967). It is not a solution of the Laplace equation in spherical coordinates and has a non-zero divergence.
4.2. The model of Maxwell (1977)
The model of Maxwell (Reference Maxwell1977) assumes an empirical form of the velocity field based on planetary cratering observations. The model assumes that the radial component 
$u_r$ is independent of 
$\theta$ and that its radial dependency is a power 
$Z$ of the radius 
$r$. Component 
$u_\theta$ is then calculated using fluid incompressibility. The velocity field is thus written as
\begin{equation} \left. \begin{aligned} & u_r=\frac{\alpha(t)}{r^Z},\\ & u_\theta=(Z-2)\frac{\sin\theta}{1+\cos\theta}\frac{\alpha(t)}{r^Z}, \end{aligned} \right\} \end{equation}
where 
$\alpha (t)$ is an arbitrary coefficient corresponding to the time-dependent flow intensity. According to Maxwell (Reference Maxwell1977) and Melosh (Reference Melosh1989), the value 
$Z=3$ gives a velocity field consistent with numerical simulations of explosion and planetary impacts.
This model allows the prediction of the experimental 
$u_{r,0}(r,t)$ (figure 10a), 
$u_{\theta,0}(r,t)$ (figure 10c) and 
$u_{\theta,1}(r,t)$ (figure 10d) using 
$Z=3$ and 
$\alpha (t)=1$. In particular, the slopes predicted by the model are very close to the experimental slopes. However, this model does not allow a degree 
$l=1$ of the radial velocity component. The main limitations of the model of Maxwell (Reference Maxwell1977) are the arbitrary choice for the model time dependency, with 
$\alpha (t)$, the fact that 
$Z$ could depend on 
$\theta$, which would yield a degree 
$l=1$ for 
$u_r$, and the fact that Maxwell's flow is not potential, which is inconsistent with the experimental results (figure 6).
4.3. The model of Leng (2001)
The model of Leng (Reference Leng2001) is similar to that of Engel (Reference Engel1967) since it uses a hemispherical crater with a radius 
$R_0(t)$ and a potential flow. The velocity potential 
$\phi$ is written as
\begin{equation} \phi={-}\frac{\dot{R}_0R_0^2}{r}, \end{equation}
which allows one to obtain the velocity components 
$u_r$ and 
$u_\theta$ of the velocity field:
\begin{equation} \left. \begin{aligned} & u_r=\frac{\dot{R}_0R_0^2}{r^2},\\ & u_\theta=0. \end{aligned} \right\} \end{equation}This velocity potential satisfies the boundary conditions and is a solution of the Laplace equation in spherical coordinates.
This model allows, in particular, the capturing of the evolution of the mean crater radius using an energy balance, although it requires one to multiply the kinetic energy and the total energy by empirical correction factors (e.g. Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). However, the velocity field has only a degree 
$l=0$ (figure 10a) on the radial component and no polar component. As for the model of Engel (Reference Engel1967), the 
$1/r^2$ slope of 
$u_{r,0}(r,t)$ is smaller than the experimental slope. The main limitations of the model of Leng (Reference Leng2001) are the hemispherical geometry and the oversimplified velocity potential which prevents a polar dependency of the radial component and a polar component of the velocity field.
4.4. The model of Bisighini et al. (2010)
The model of Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) assumes an expanding spherical crater able to translate vertically over time, with a radius 
$R_0(t)$ and a vertical position of the crater barycentre 
$z_c(t)$. This allows the definition of a velocity potential 
$\phi$ which corresponds to the superposition between the radial expansion of the crater and the flow past a translating sphere. This potential satisfies the boundary conditions and the Laplace equation in spherical coordinates. In the moving sphere coordinate system 
$(r',\theta ')$, it is written
\begin{equation} \phi={-}\frac{\dot{R}_0R_0^2}{r'}-\dot{z}_c r'\left(1-\frac{R_0^3}{2r'^3}\right)\cos\theta', \end{equation}
with components 
$u_r$ and 
$u_\theta$ of the velocity field:
\begin{equation} \left. \begin{aligned} u_r & =\frac{\dot{R}_0R_0^2}{r'^2}-\left(1-\frac{R_0^3}{r'^3}\right)\dot{z}_c\cos\theta',\\ u_\theta & =\left(1+\frac{R_0^3}{2r'^3}\right)\dot{z}_c\sin\theta' . \end{aligned} \right\} \end{equation}
Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) then use an unsteady Bernoulli equation to determine the evolution of the sphere radius and position over time. To compare the model of Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) with our experimental data, we need to calculate the corresponding velocity field in the fixed frame of reference by adding the velocity of the crater barycentre 
$\dot {z}_c (\cos \theta, -\sin \theta )$ to (4.7), and expressing 
$r'$ and 
$\theta '$ as functions of 
$r$ and 
$\theta$ (
$r'=\sqrt {r^2+z_c^2-2 z_c r \cos \theta }$, 
$\cos \theta '=(r\cos \theta - z_c)/r'$, 
$\sin \theta '=r\sin \theta /r'$).
The velocity field has an 
$l=0$ (figure 10a) and an 
$l=1$ (figure 10b) radial component, as well as an 
$l=0$ (figure 10c) and an 
$l=1$ (figure 10d) polar component. The coefficients are calculated using 
$z_c=0$ and 
$\dot {z}_c=0.2U_i$, which correspond to typical values during crater opening (e.g. figure 5). This model explains relatively well the shape of the crater (e.g. Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010, figure 17), and the key tendencies of the experimental components of the velocity field. However, the model of Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) strongly constrains the geometry of the crater, as well as the related velocity potential definition. As in the models of Engel (Reference Engel1967) and Leng (Reference Leng2001), the 
$1/r^2$ slope of 
$u_{r,0}$ is smaller than the experimental slope.
4.5. Towards a new model
In all models, either the geometry of the velocity field (Engel Reference Engel1967; Maxwell Reference Maxwell1977; Leng Reference Leng2001) or the shape of the cavity (Engel Reference Engel1967; Leng Reference Leng2001; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010) is imposed. This leads in particular to an incorrect radial dependency of 
$u_r$, with an exponent much larger in the experiments than in the models, except for the model of Maxwell (Reference Maxwell1977) model where the radial dependency is arbitrarily imposed by the parameter 
$Z$. The experimental observation that the radial velocity field decreases with 
$r$ faster than 
$1/r^2$ is unexpected since it suggests that the flow component associated with an isotropic expansion of the cavity (
$\propto 1/r^2$) is not dominant. New models are thus required to explain the geometry of the experimental velocity field, as well as the evolution of the non-hemispherical shape of the cavity. In the following section, we develop a semi-analytical model based on a Legendre polynomial expansion of an unsteady Bernoulli equation, coupled with a kinematic boundary condition at the crater boundary.
5. Legendre polynomial model
In this model, we assume that the fluid is inviscid (i.e. 
$\mu =0$), incompressible (i.e. 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}=0$) and that the flow is irrotational (i.e. 
$\boldsymbol {\nabla }\times \boldsymbol {u}=0$). This means that the flow is potential and satisfies the Laplace equation 
$\nabla ^2\phi =0$, where 
$\phi$ is the velocity potential defined as 
$\boldsymbol {u}=\boldsymbol {\nabla }\phi$. In the spherical coordinate system 
$(r,\theta,\varphi )$, assuming an axisymmetric flow, the solution of the Laplace equation is written as
\begin{equation} \phi(r,\theta,t)=\sum_{n=0}^{+\infty}\frac{\phi_n(t)}{r^{n+1}}P_n(\cos\theta), \end{equation}
where 
$\phi _n(t)$ are time-dependent coefficients and 
$P_n(x)$ are the standard Legendre polynomials, orthogonal on 
$[-1,1]$. The components 
$u_r$ and 
$u_\theta$ of the velocity field are then
\begin{equation} \left. \begin{aligned} u_r(r,\theta,t) & =\frac{\partial \phi}{\partial r}= \sum_{n=0}^{+\infty}-(n+1)\frac{\phi_n(t)}{r^{n+2}}P_n(\cos\theta), \\ u_\theta(r,\theta,t) & =\frac{1}{r}\frac{\partial \phi}{\partial\theta}= \sum_{n=0}^{+\infty}\frac{\phi_n(t)}{r^{n+2}}\frac{\partial P_n(\cos\theta)}{\partial\theta}. \end{aligned} \right\} \end{equation}We also assume a non-hemispherical crater, where the shape of the cavity is decomposed on a set of shifted Legendre polynomials (2.2).
Since we assume that the fluid is inviscid and a potential flow, the flow is governed by an unsteady Bernoulli equation:
\begin{equation} \frac{\partial\phi}{\partial t} + \frac{1}{2}u^2 - gz + \frac{p}{\rho} = \text{constant}, \end{equation}
where 
$\rho$ is the fluid density, 
$u$ is the norm of the velocity, 
$p$ is the pressure, 
$g$ is the acceleration due to gravity and 
$z$ is the vertical coordinate below the initial fluid surface. This equation is constant in the entire fluid domain. Far from the crater, 
$u \to 0$, 
$\phi \to 0$ and the pressure is hydrostatic 
$p(z)=p_0+\rho g z$, where 
$p_0$ is the atmospheric pressure. This means that the constant is equal to 
$p_0/\rho$.
At the crater boundary, i.e. at 
$r=R(\theta,t)$ (equation (2.2)), the Young–Laplace equation is
\begin{equation} p(R)-p_0=\sigma C, \end{equation}
where 
$C(\theta,t)$ is the mean local curvature of the interface and 
$\sigma$ the surface tension. In cylindrical coordinates, the curvature is
\begin{equation} C(\theta,t)=\frac{R^2+2\left(\dfrac{\partial R}{\partial\theta}\right)^2-R\dfrac{\partial^2R}{\partial\theta^2}}{\left[R^2+\left(\dfrac{\partial R}{\partial\theta}\right)^2\right]^{3/2}}. \end{equation}The Bernoulli equation at the crater boundary is thus written as
\begin{equation} \left(\dfrac{\partial\phi}{\partial t}\right)_{r=R}+\dfrac{1}{2}u(R)^2-gR\cos\theta+\dfrac{\sigma}{\rho}C=0. \end{equation}We also use a kinematic boundary condition at the crater boundary:
\begin{equation} \dfrac{\partial R}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} R=\dfrac{\partial R}{\partial t}+u_\theta(R)\dfrac{1}{R}\dfrac{\partial R}{\partial \theta}=u_r(R). \end{equation} Equations (5.6) and (5.7) are made dimensionless using the scaling laws for the crater opening time scale 
$\tilde {t}_{max}$ (3.1) and the maximum crater radius 
$\tilde {R}_{max}$ (3.2), which gives the following partial differential equation system:
\begin{align}
\left. \begin{aligned} \left(\frac{\partial\phi^*}{\partial
t^*}\right)_{r^*=R^*} &
={-}\frac{1}{2}u^*(R^*)^2+\frac{1}{4}\mathrm{B}\left(\!\frac{1}{2},\frac{5}{8}\!\right)^2
\xi R^*\cos\theta-\frac{1}{8}\sqrt{\frac{3}{2}}\mathrm{B}\left(\!\frac{1}{2},\frac{5}{8}\!\right)^2
\frac{\xi}{\sqrt{\varPhi}}\frac{\sqrt{Fr}}{We}C^*,\\
\frac{\partial R^*}{\partial t^*} &
=u^*_r(R^*)-u^*_\theta(R^*)\frac{1}{R^*}\frac{\partial
R^*}{\partial \theta}, \end{aligned} \right\}
\end{align}
where the star notation denotes quantities made dimensionless with 
$\tilde {R}_{max}$ and 
$\tilde {t}_{max}$, e.g. 
$t^*=\tilde {t}/\tilde {t}_{max}$.
We solve this differential equation system (5.8) by expanding the velocity potential (5.1) up to degree 
$n_{max}=2$:
\begin{equation} \phi^*(r^*,\theta,t^*)=\frac{\phi^*_0(t^*)}{{r^*}}+\frac{\phi^*_1(t^*)\cos\theta}{{r^*}^2}+\frac{\phi^*_2(t^*)(3\cos^2\theta-1)}{2{r^*}^3}. \end{equation}The components of the velocity field are then written (equation (5.2))
\begin{equation} \left. \begin{aligned} u^*_r(r^*,\theta,t^*) & ={-}\frac{\phi^*_0(t^*)}{{r^*}^2}-\frac{2 \phi^*_1(t^*)\cos\theta}{{r^*}^3}-\frac{3\phi^*_2(t^*)(3\cos^2\theta-1)}{2{r^*}^4},\\ u^*_\theta(r^*,\theta,t^*) & ={-}\frac{\phi^*_1(t^*)\sin\theta}{{r^*}^3}-\frac{3\phi^*_2(t^*)\sin\theta\cos\theta}{{r^*}^4}. \end{aligned} \right\} \end{equation}
We also expand the crater boundary position (2.2) up to degree 
$k_{max}=1$:
\begin{equation} R^*(\theta,t^*)=R^*_0(t^*)+R^*_1(t^*)(2\cos\theta-1). \end{equation}
Note that the crater position 
$R^*(\theta,t^*)$ is written as a sum of shifted Legendre polynomials, while the velocity potential 
$\phi ^*(r^*,\theta,t^*)$ is a sum of standard Legendre polynomials.
We then project the differential equation system (5.8) on a set of shifted Legendre polynomials 
$\bar {P}_m$ up to degree 
$m_{max}=2$ for the Bernoulli equation and degree 
$m_{max}=1$ for the kinematic boundary condition. The projection of a function 
$X$ is written as
\begin{equation} \langle X,\bar{P}_m \rangle = (2m+1)\int_{0}^{{\rm \pi}/2}X\bar{P}_m(\cos\theta)\sin\theta\,\mathrm{d}\theta. \end{equation}
We simplify the equations by expanding the Bernoulli equation and the kinematic boundary condition to the third and the fourth order in 
$R^*_1$. We obtain a system of five equations with five unknowns: 
$\phi ^*_0(t^*)$, 
$\phi ^*_1(t^*)$, 
$\phi ^*_2(t^*)$, 
$R^*_0(t^*)$ and 
$R^*_1(t^*)$ (equations (A1)–(A5)).
The general equation system (5.8) and its projection (A1)–(A5) may be further simplified. The third term on the right-hand side of the Bernoulli equation (in (5.8)) corresponds to surface tension effects associated with the curvature of the air–water interface. If this term is neglected, which corresponds to 
$\sqrt {Fr}/We \ll 1$, (5.8) then simplifies as
\begin{equation} \left. \begin{aligned} \left(\frac{\partial\phi^*}{\partial t^*}\right)_{r^*=R^*} & ={-}\frac{1}{2}u^*(R^*)^2+\frac{1}{4}\mathrm{B}\left(\frac{1}{2},\frac{5}{8}\right)^2 \xi R^*\cos\theta,\\ \frac{\partial R^*}{\partial t^*} & =u^*_r(R^*)-u^*_\theta(R^*)\frac{1}{R^*}\frac{\partial R^*}{\partial \theta}. \end{aligned} \right\} \end{equation}
In our experiments, 
$\sqrt {Fr}/We$ is two to three times larger for case A (
$1.0 \times 10^{-1}$) than for cases B, C and D (
$4.9 \times 10^{-2}$, 
$3.9 \times 10^{-2}$ and 
$3.3 \times 10^{-2}$, respectively). This is consistent with the surface tension argument used to explain the difference between case A and the other cases (§ 3). Since 
$\xi$ is independent of 
$Fr$ and 
$We$ in our experimental range (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022), this normalised equation system without surface tension is independent of the impact parameters and may be used to provide a predictive model.
The general and the simplified equation systems are solved numerically as initial value problems, using a differential equation solver. The solution thus depends on the choice of initial conditions. On the one hand, we can solve the equation systems separately for each experiment. The initial conditions are defined at 
$\tilde {t}=1$, which corresponds to an advection time of the impacting drop, and fitted on each experiment by using a joint least-squares inversion of the five experimental coefficients over the entire time series. On the other hand, we can solve the equation systems at the same time for all the experiments. The initial conditions are also defined at 
$\tilde {t}=1$ but fitted simultaneously on all the experiments using the joint least-squares inversion over the entire time series. This method allows the definition of a unique set of initial conditions that may be used in a predictive model. In both cases, the fitting procedure is motivated by the sensitivity of the model to the initial conditions used. A slight modification of the initial conditions may change significantly the time evolution of the coefficients. This sensitivity might be related to the exact impact conditions, including a possible variability in the contact dynamics with the surface of the pool and in the shape of the drop upon impact. Furthermore, the sensitivity to initial conditions might be amplified by the truncation of the crater shape and of the velocity potential expansion, which is probably insufficient to model properly the early evolution of the crater. This sensitivity is investigated in more detail in Appendix B.
We now define two models using different systems of equations and definitions of initial conditions. The first model, referred to as the general model, accounts for surface energy effects and uses the general equation system (5.8) and initial conditions fitted on single experiments. This means that the number of sets of initial conditions is equal to the number of experiments. For example, the initial conditions of a given experiment in case B are 
$\phi ^*_0(1)=-0.07 \pm 0.02$, 
$\phi ^*_1(1)=-0.07 \pm 0.02$, 
$\phi ^*_2(1)=0.009 \pm 0.003$, 
$R^*_0(1)=0.41 \pm 0.03$ and 
$R^*_1(1)=-0.28 \pm 0.02$. Uncertainties on the coefficients correspond to 
$1-\sigma$ standard deviations on the parameters in the least-squares inversion. The initial conditions of all the experiments are presented in Appendix B. The second model, referred to as the simplified model, uses the simplified equation system, without surface tension and independent of the impact parameters (5.13), as well as initial conditions fitted on all the experiments. The reference set of initial conditions is
\begin{align} \left. \begin{aligned}
& \phi^*_0(1)={-}0.21 \pm 0.01, & \phi^*_1(1)=0.002 \pm
0.005, & \quad\phi^*_2(1)=0.0004 \pm 0.0005,\\ & R^*_0(1)=0.29
\pm 0.02, & R^*_1(1)={-}0.39 \pm 0.02. \end{aligned}
\right\} \end{align}
Given the uncertainties, this set of initial conditions can be further simplified by using 
$\phi ^*_1(1)=\phi ^*_2(1)=0$, which corresponds to an initial velocity field given by 
$(u^*_r=-\phi ^*_0(1)/{r^*}^2, u^*_\theta =0)$. The physical interpretation of these initial conditions should be investigated in the future. It probably involves the contact dynamics between the drop and the pool and the early evolution of the crater. Nonetheless, the simplified model is a predictive model, independent of the impact parameters, that can be used to predict the crater and velocity field evolution within the range of 
$Fr$ and 
$We$ covered by our experiments. However, we anticipate the model to show predictability limitations outside of this range, in particular at low 
$Fr$ and 
$We$ in the bubble entrapment region (e.g. Pumphrey & Elmore Reference Pumphrey and Elmore1990), due to the neglected surface tension term and more generally to the relatively low degree of truncation used in our model.
Figure 11 compares the experimental coefficients 
$\phi ^*_0$ (figure 11a), 
$\phi ^*_1$ (figure 11b) and 
$\phi ^*_2$ (figure 11c) of the velocity potential and the experimental coefficients 
$R^*_0$ (figure 11d) and 
$R^*_1$ (figure 11e) of the crater shape with the coefficients obtained with the general (coloured solid lines) and the simplified (black solid lines) models. We determine the experimental velocity potential coefficients from the experimental velocity field using a joint least-squares inversion of the radial and the polar components (5.2). We also obtain the experimental crater shape coefficients by fitting the crater boundary position with the shifted Legendre polynomial expansion (2.2), using the method described in § 2.4.1.
Figure 11. Coefficients 
$\phi ^*_0$ (a), 
$\phi ^*_1$ (b), 
$\phi ^*_2$ (c), 
$R^*_0$ (d) and 
$R^*_1$ (e) as a function of time normalised by the opening time scale of the crater 
$\tilde {t}/\tilde {t}_{max}$, in the four cases. For each case, the different types of markers correspond to different experiments. The coloured solid lines give the solution of the general model, where the initial conditions are fitted on a single experiment of the corresponding case. The black solid lines give a solution to the simplified model, where the initial conditions are fitted simultaneously on all the experiments. The black dashed lines give a solution to the simplified model, where the initial conditions are modified by 
$\pm 25\,\%$ with respect to their reference value.
The models capture well the evolution of the velocity potential (figure 11a–c) for all cases. In detail, the models are dominated by 
$\phi ^*_1$ and are slightly less accurate when it comes to fit 
$\phi ^*_2$, as expected since this corresponds to velocity fluctuations on smaller scales. These results are consistent with the good agreement between the simplified model and the experimental velocity coefficients of figure 10, in particular regarding the slope of 
$u_{r,0}(r,t)$. Although 
$u_{r,0}(r,t)$ remains less steep than in the experiments, it decreases significantly faster than 
$1/r^2$. The models also capture well the evolution of the crater shape (figure 11d,e). Note that 
$\tilde {R}_{max}$ slightly overestimates the experimental maximum crater radius, with maximum 
$R^*_0$ systematically smaller than 1. This can be explained by the neglected surface energy in the energy balance (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022). In detail, 
$R^*_1$ is slightly underestimated and changes at a higher rate than experimental data when 
$\tilde {t}/\tilde {t}_{max} \lesssim 0.4$ and 
$\tilde {t}/\tilde {t}_{max} \gtrsim 1.7$. This corresponds respectively to the early opening of the crater and the end of crater collapse, including the formation of the central jet, where an expansion of 
$R$ to a higher degree (at least 
$k=2$) would be required to model the observed degree of deformation of the cavity (e.g. figure 4c).
Note that the predictive model, using the simplified equation system (5.13) with the reference set of initial conditions (5.14), is particularly in good agreement with the experimental data. The sensitivity of the simplified model to the initial conditions is illustrated with two solutions where the initial conditions have been modified by 
$\pm 25\,\%$ with respect to their reference value (figure 11, black dashed lines).
Although case A is slightly different from cases B, C and D due to surface tension effects (see § 3), the models capture properly the general cratering dynamics. In detail, 
$\phi ^*_2$ and 
$R^*_1$ are significantly underestimated when 
$0.5 \lesssim \tilde {t}/\tilde {t}_{max} \lesssim 1.4$, as expected since the models do not account for the capillary wave propagation responsible for this cavity deformation.
Figure 12 compares snapshots of the radial (figure 12a,b) and polar (figure 12c,d) components of the experimental velocity field (figure 12a,c) with the components calculated from the predictive simplified model (figure 12b,d), in case B. The comparison is conducted at different times during the opening stage (i), just before the crater reaches its maximum size (ii), and during the closing stage (iii). This illustrates that the velocity fields from the simplified model and the experiment are very similar during all stages of the cratering process. The differences observed are mainly in the magnitude of the velocity, in particular close to the crater and the initial water surface (
$\theta = \pm {\rm \pi}/2$). Similar results are obtained in the other cases. The good agreement between the experimental velocity field and the simplified model shows that the truncation used in the model (degree 
$k=1$ in shifted Legendre polynomials for 
$R^*$ and degree 
$n=2$ in Legendre polynomials for 
$\phi ^*$) is sufficient to accurately capture the flow dynamics.
Figure 12. Radial (a,b) and polar (c,d) components of the velocity field from experimental data (a,c) and from the simplified model (b,d), in case B. The snapshots correspond to times when the crater is opening (i), when the crater is almost at its maximum size (ii) and when the crater is closing (iii). The solid green lines correspond to the experimental crater boundary.
Figure 13 compares the crater shape obtained in a backlight experiment similar to case B (
$Fr=442$) with the crater boundary position calculated from the predictive simplified model. The crater shape is well captured by the model, consistent with the results of figure 11(d,e). In detail, at the very beginning of the crater opening stage (figure 13i), the model overestimates the width of the crater and does not capture accurately the flat-bottomed shape of the cavity. During the crater opening stage and the beginning of the crater collapse stage (figure 13iii–v), the model slightly underestimates the crater depth and width, consistent with the coefficients of figure 11(d,e). Finally, when the crater collapses (figure 13vi), the model shows the central jet initiation, although it visibly lacks higher degrees to account for the vertical walls of the cavity. Figure 13 also compares the experimental velocity field obtained in case B with the velocity field obtained from the simplified model. The comparison shows a good agreement between the two, which is consistent with the analysis of figure 12.
Figure 13. Time evolution of the crater shape obtained in a backlight experiment similar to case B (
$Fr=442$). The solid green lines correspond to the crater boundary obtained from the simplified model. The black arrows correspond to the experimental velocity field, normalised by its maximum value in each snapshot. The grey arrows correspond to the velocity field obtained from the simplified model, normalised by its maximum value in each snapshot.
6. Conclusion
In this paper, we analyse quantitatively the velocity field around the crater produced by the impact of a liquid drop onto a deep liquid pool. Using new high-resolution PIV measurements, we obtain simultaneously the evolution of the velocity field around the cavity and the crater shape. We found that the shape of the cavity and the velocity field seem to be independent of 
$Fr$ and 
$We$ at a given 
$t/t_{max}$, when these two dimensionless numbers are large enough (cases B, C and D). The velocity field is dominated by the degrees 0 and 1 in terms of shifted Legendre polynomials, with the degree 0 of the radial component 
$u_{r,0}(r,t)$ decreasing faster than 
$1/r^2$. Furthermore, the radial component of the velocity field is dominated by the degree 1 in terms of standard Legendre polynomials. This is not inconsistent with the growth of the crater because the degree 1 of the radial component has a non-zero average over a hemisphere. The experiments also show significant contributions from the degree 2, in particular when the crater is strongly deformed. This is possibly related to the non-hemispherical shape (degree 1) of the cavity. We also found that the velocity field does not vanish when the crater reaches its maximum size.
In previous velocity models (Engel Reference Engel1967; Maxwell Reference Maxwell1977; Leng Reference Leng2001; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010), strong constraints were imposed on the crater shape and/or on the velocity field. They were unable to explain the properties observed in our experimental measurements, in particular the radial dependency of the radial component of the velocity field and the evolution of the shape of the cavity. We thus developed a semi-analytical model based on a Legendre polynomial expansion of an unsteady Bernoulli equation, coupled with a kinematic boundary condition at the crater boundary. Assuming that the surface tension term involved in the Bernoulli equation is negligible, we define a simplified model, independent of the impact parameters, that can predict the evolution of the crater shape and of the velocity field within the range of 
$Fr$ and 
$We$ numbers covered in our experiments. Although the model is sensitive to the initial conditions, it remains predictive by using a unique set of fitted initial conditions. In particular, the model can capture the initiation of the central jet. However, one intrinsic limitation of the model is that it assumes the cavity radius to be a bijective function of 
$\theta$. While this assumption is true during the opening stage and part of the crater closing stage, including the central jet initiation, it eventually fails when the central jet reaches a critical height, since the air–water interface can be crossed twice at a given 
$\theta$. The model can therefore not be used to describe the full growth of the central jet.
Acknowledgements
We thank M. Moulin for his help with the design and construction of the experimental apparatus. We thank three anonymous reviewers for their valuable comments which significantly improved the manuscript.
Funding
This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant number 716429). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). Partial funding for this research was provided by the Center for Matter at Atomic Pressure (CMAP), a National Science Foundation (NSF) Physics Frontier Center, under award PHY-2020249. Any opinions, findings, conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of the National Science Foundation.
Declaration of interests
The authors report no conflict of interest.
Author contributions
V.L. and R.D. designed the study, derived the model and contributed to analysing data and reaching conclusions. V.L. conducted the experiments. V.L. and R.D wrote the paper.
Appendix A. Equations of the Legendre polynomial model
The Legendre polynomial model equations correspond to the projection of (5.8) up to degree 
$m_{max}=1$ for the kinematic boundary condition and up to degree 
$m_{max}=2$ for the Bernoulli equation. The projected boundary conditions and Bernoulli equations are respectively expanded to the fourth and the third order in 
$R^*_1$. The boundary condition is then written as
\begin{align}
\dot{R}^*_0&=\frac{1}{6 {R^*_0}^8}\left\{{-}6 {\phi^*_0}
{R^*_0}^2 \left({R^*_0}^2
{R^*_1}^2+{R^*_0}^4+{R^*_1}^4\right)-2 {\phi^*_1} {R^*_0}
\left({R^*_0}^3 {R^*_1}+10 {R^*_0}^2 {R^*_1}^2+6 {R^*_0}
{R^*_1}^3\right.\right.\nonumber\\ &\quad \left.\left.+3
{R^*_0}^4+21 {R^*_1}^4\right)+3 {\phi^*_2} {R^*_1} (3
{R^*_0}-{R^*_1}) \left({R^*_0}^2+3
{R^*_1}^2\right)\right\}\nonumber\\ &\quad
+O({R^*_1}^5),
\end{align}
\begin{align}\dot{R}^*_1&=\frac{1}{140
{R^*_0}^8}\left\{14 {R^*_0}^3 {R^*_1} [20 {R^*_0}
({\phi^*_0} {R^*_0}+2 {\phi^*_1})+9 {\phi^*_2}]\right.\nonumber\\
&\quad +6 {R^*_0}
{R^*_1}^3 [56 {R^*_0} ({\phi^*_0} {R^*_0}+5 {\phi^*_1})+75
{\phi^*_2}] +84 {R^*_0}^2
{R^*_1}^2 (2 {\phi^*_1} {R^*_0}\nonumber\\
&\left.\quad -15 {\phi^*_2})+15 {R^*_1}^4
(68 {\phi^*_1} {R^*_0}-189 {\phi^*_2})-35 {R^*_0}^4 (4
{\phi^*_1} {R^*_0}+9 {\phi^*_2})\right\} +O({R^*_1}^5),
\end{align}and the Bernoulli equation is written as
\begin{align}
0&=\frac{1}{1680 {R^*_0}^{11}}\left\{28 {R^*_0}^2
\left[{R^*_0}^3 \left(20 {R^*_0}^3 \dot{\phi}^*_0 \left(3
{R^*_0}^2+{R^*_1}^2\right)+10 {\phi^*_0}^2 \left(3
{R^*_0}^2+10
{R^*_1}^2\right)\right.\right.\right.\nonumber\\
&\quad \left.\left.\left.+2 {R^*_0} \dot{\phi}^*_1 \left({-}10
{R^*_0}^2 {R^*_1}+15 {R^*_0} {R^*_1}^2+15 {R^*_0}^3-12
{R^*_1}^3\right)-3 {R^*_1} \dot{\phi}^*_2 \left({-}4
{R^*_0} {R^*_1}+15
{R^*_0}^2\right.\right.\right.\right.\nonumber\\
&\quad \left.\left.\left.\left.+30 {R^*_1}^2\right)\right)+20
{\phi^*_0} {\phi^*_1} {R^*_0} \left({-}5 {R^*_0}^2
{R^*_1}+15 {R^*_0} {R^*_1}^2+3 {R^*_0}^3-21
{R^*_1}^3\right)\right.\right.\nonumber\\
&\quad \left.\left.\left.+6 {\phi^*_1}^2 \left({-}15 {R^*_0}^2
{R^*_1}+77 {R^*_0} {R^*_1}^2+10 {R^*_0}^3-84
{R^*_1}^3\right)\right]+168 {\phi^*_2} {R^*_0} \left[{-}9
{R^*_0}^2 {R^*_1} (5 {\phi^*_0} {R^*_0}\right.\right.\right.\nonumber\\
&\quad \left.\left.+7 {\phi^*_1})+7
{R^*_0} {R^*_1}^2 (3 {\phi^*_0}
{R^*_0}+28 {\phi^*_1})-36 {R^*_1}^3 (7 {\phi^*_0} {R^*_0}+13
{\phi^*_1})+15 {\phi^*_1} {R^*_0}^3\right]\right.\nonumber\\
&\quad \left.+72 {\phi^*_2}^2
\left[{-}70 {R^*_0}^2 {R^*_1}+558 {R^*_0}
{R^*_1}^2+35
{R^*_0}^3-720 {R^*_1}^3\right]\right.\nonumber\\
&\quad \left.\left.-\frac{35}{\sqrt{\varPhi}
We}\mathrm{B}\left(1/2,5/8\right)^2 \xi {R^*_0}^6
\left[{-}3 \sqrt{6} \sqrt{Fr} {R^*_0}^2
{R^*_1}^2 +3 \sqrt{6} \sqrt{Fr} {R^*_0} {R^*_1}
\left({R^*_1}^2-{R^*_0}^2\right)\right.\right.\right.\nonumber\\
&\quad \left.\left.-3 \sqrt{6} \sqrt{Fr}
{R^*_0}^4+2 \sqrt{\varPhi} We {R^*_0}^5 {R^*_1}+6
\sqrt{\varPhi} We {R^*_0}^6\right]\right\}+O({R^*_1}^4),
\end{align}
\begin{align}
0&=\frac{1}{280 {R^*_0}^{11}}\left\{{-}14 {R^*_0}^2 \left[4
{R^*_0}^2 {R^*_1} \left({R^*_0}^3 \dot{\phi}^*_0 \left(5
{R^*_0}^2+3 {R^*_1}^2\right)+10 {\phi^*_0}^2
\left({R^*_0}^2+3
{R^*_1}^2\right)\right)\right.\right.\nonumber\\
&\quad\left.\left.-20 {\phi^*_0} {\phi^*_1} {R^*_0}
\left({-}5 {R^*_0}^2 {R^*_1}+9 {R^*_0}
{R^*_1}^2+{R^*_0}^3-21 {R^*_1}^3\right)+2 {R^*_0}^4
\dot{\phi}^*_1 \left(10 {R^*_0}^2 {R^*_1}-9 {R^*_0}
{R^*_1}^2\right.\right.\right. \nonumber\\
&\quad\left.\left.\left.-5 {R^*_0}^3+12 {R^*_1}^3\right)+3
{\phi^*_1}^2 \left(44 {R^*_0}^2 {R^*_1}-63 {R^*_0}
{R^*_1}^2-5 {R^*_0}^3+256
{R^*_1}^3\right)\right]\right.\nonumber\\ &\quad\left.+42
{\phi^*_2} {R^*_0} \left[3 {\phi^*_0} {R^*_0} \left({-}4
{R^*_0}^2 {R^*_1}+63 {R^*_0} {R^*_1}^2+5 {R^*_0}^3-32
{R^*_1}^3\right)+2 {\phi^*_1} \left({-}49 {R^*_0}^2
{R^*_1}\right.\right.\right.\nonumber\\
&\quad\left.\left.\left.+156 {R^*_0} {R^*_1}^2+9
{R^*_0}^3-396 {R^*_1}^3\right)\right]+6 {R^*_0}^5
\dot{\phi}^*_2 \left[{-}14 {R^*_0}^2 {R^*_1}+126 {R^*_0}
{R^*_1}^2\right.\right.\nonumber\\
&\quad\left.\left.+35 {R^*_0}^3-40
{R^*_1}^3\right]+9 {\phi^*_2}^2 \left[{-}496 {R^*_0}^2 {R^*_1}+864 {R^*_0}
{R^*_1}^2+35 {R^*_0}^3-4960
{R^*_1}^3\right]\right.\nonumber\\
&\quad\left.-\frac{1120 {\rm \pi}}{\sqrt{\varPhi} We }\frac{\varGamma
\left(5/8\right)^2}{\varGamma \left(1/8\right)^2} \xi
{R^*_0}^6 ({R^*_0}+{R^*_1}) \left[3 \sqrt{6} \sqrt{Fr}
{R^*_0} {R^*_1}^2+2 \sqrt{\varPhi} We
{R^*_0}^5\right]\right\}\nonumber\\
&\quad+O({R^*_1}^4), \end{align}
\begin{align} 0&=\frac{1}{336 {R^*_0}^{11}}\left\{{R^*_0}^2 \left[{R^*_0}^3 \left(4 \left(56 {R^*_1}^2 \left({R^*_0}^3 \dot{\phi}^*_0+5 {\phi^*_0}^2\right)-4 {R^*_0} {R^*_1} \dot{\phi}^*_1 \left({-}21 {R^*_0} {R^*_1}+14 {R^*_0}^2\right.\right.\right.\right.\right.\nonumber\\ &\quad\left.\left.\left.\left.\left.+24 {R^*_1}^2\right)+3 \dot{\phi}^*_2 \left({-}42 {R^*_0}^2 {R^*_1}+22 {R^*_0} {R^*_1}^2+7 {R^*_0}^3-120 {R^*_1}^3\right)\right)\right.\right.\right.\nonumber\\ &\quad\left.\left.\left.-\frac{7}{\sqrt{\varPhi} We}\mathrm{B}\left(1/2,5/8\right)^2 \xi {R^*_1} \left(3 \sqrt{6} \sqrt{Fr} {R^*_0}^3 {R^*_1}-21 \sqrt{6} \sqrt{Fr} {R^*_0}^2 {R^*_1}^2\right.\right.\right.\right.\nonumber\\ &\quad\left.\left.\left.\left.+4 \sqrt{\varPhi} We {R^*_0}^6\right)\right)-1120 {\phi^*_0} {\phi^*_1} {R^*_0} {R^*_1} \left({-}3 {R^*_0} {R^*_1}+{R^*_0}^2+6 {R^*_1}^2\right)+84 {\phi^*_1}^2 \left({-}12 {R^*_0}^2 {R^*_1}\right.\right.\right.\nonumber\\ &\quad\left.\left.\left.+67 {R^*_0} {R^*_1}^2+{R^*_0}^3-96 {R^*_1}^3\right)\right]+84 {\phi^*_2} {R^*_0} \left[3 {\phi^*_0} {R^*_0} \left({-}12 {R^*_0}^2 {R^*_1}+11 {R^*_0} {R^*_1}^2+{R^*_0}^3\right.\right.\right.\nonumber\\ &\quad\left.\left.-96 {R^*_1}^3\right)+2 {\phi^*_1} ({R^*_0}-6 {R^*_1}) \left({-}9 {R^*_0} {R^*_1}+3 {R^*_0}^2+46 {R^*_1}^2\right)\right]+18 {\phi^*_2}^2 \left[{-}148 {R^*_0}^2 {R^*_1}\right.\nonumber\\ &\quad\left.\left.+1116 {R^*_0} {R^*_1}^2+23 {R^*_0}^3-1860 {R^*_1}^3\right]\right\}\nonumber\\ &\quad+O({R^*_1}^4). \end{align}
The simplified version of the equation system (5.13) can be obtained by using 
$\sqrt {Fr}/We=0$ in (A1)–(A5).
Appendix B. Initial conditions of the Legendre polynomial model
Figure 14 shows the initial conditions of the general model, obtained by fitting individually the experiments, and of the simplified model, obtained by fitting all the experiments simultaneously. They are both defined at 
$\tilde {t}=1$ and use a joint least-squares inversion of the experimental coefficients over the entire time series. Uncertainties on the initial conditions correspond to 
$1-\sigma$ standard deviations on the parameters in the least-squares inversion.
Figure 14. Initial conditions (at 
$\tilde {t}=1$) of the Legendre polynomial models, for the coefficients 
$\phi ^*_0$ (a), 
$\phi ^*_1$ (b), 
$\phi ^*_2$ (c), 
$R^*_0$ (d) and 
$R^*_1$ (e), obtained using a joint least-squares inversion of the experimental coefficients over the entire time series. The coloured markers correspond to the general model, where initial conditions are obtained by fitting experiments individually. The black circles correspond to the simplified model, where initial conditions are obtained by fitting all the experiments simultaneously. Uncertainties correspond to 
$1-\sigma$ standard deviations on the parameters in the least-squares inversion.
At low 
$\sqrt {Fr}/We$, corresponding to high 
$Fr$ and 
$We$ numbers (cases B, C, D), the dispersion of the initial conditions is larger than the uncertainties associated with the least-squares inversion, whereas at higher 
$\sqrt {Fr}/We$, corresponding to moderate 
$Fr$ and 
$We$ numbers (case A), the initial conditions are clustered within the inversion uncertainties. This dispersion at higher 
$Fr$ and 
$We$ suggests a higher variability of the crater shape and of the velocity field upon impact. This might be related to a greater sensitivity to the exact impact conditions, possibly including variability in the contact dynamics with the surface of the pool and in the shape of the drop upon impact. Furthermore, we do not find any secondary dependency on 
$Fr$ or 
$We$. Finally, the initial conditions of the simplified model, obtained by fitting all the experiments simultaneously, are similar to the initial conditions obtained by fitting individually the experiments.
The relatively large dispersion observed for a given case (except for case A) indicates that the model is sensitive to the initial conditions. For example, a change in all the initial conditions by 
$\pm 25\,\%$ gives a significantly modified evolution of the coefficients over time (figure 11, black dashed lines). In order to further investigate this initial condition sensitivity, we conducted a quantitative test on the simplified model. Figure 15 shows the relative change of the model coefficients with respect to the simplified model, as a result of an individual modification of a single initial condition from the reference value defined in (5.14). The relative change 
$\delta X$ is defined as the absolute change in 
$X=\{\phi ^*_0, \phi ^*_1, \phi ^*_2, R^*_0, R^*_1\}$, 
$X-X_{ref}$, scaled by the root mean square of the simplified model 
$\mathrm {r.m.s.}(X_{ref})$. We choose to scale the absolute change by the root mean square of the simplified model to ensure a non-diverging value of the relative change when 
$X_{ref} \to 0$. Note that this sensitivity test only investigates the role of independent parameter modifications. Coupled modifications of the initial conditions (as in figure 11, black dashed lines) might amplify significantly the changes in the evolution of the coefficients.
Figure 15. Sensitivity of the simplified model to individual changes in the initial conditions. (a–e) Each row corresponds to the relative change 
$\delta X(1)=[X(1)-X_{ref}(1)]/X_{ref}(1)$ (in 
$\%$) of a single initial condition 
$X(1)=\{\phi ^*_0(1), \phi ^*_1(1), \phi ^*_2(1), R^*_0(1), R^*_1(1)\}$ with respect to the reference value 
$X_{ref}(1)$, defined in (5.14). Each column corresponds to the relative change 
$\delta X=(X-X_{ref})/\mathrm {r.m.s.}(X_{ref})$ (in 
$\%$) of the coefficient 
$X$ over the entire time series, for a given initial condition modification. The relative change 
$\delta X$ is defined as the absolute change in 
$X$, 
$X-X_{ref}$, scaled by the root mean square of the simplified model 
$\mathrm {r.m.s.}(X_{ref})$.
Within the range of parameter modifications (by 
$\pm 40\,\%$), the coefficients are generally more influenced by modifications of the initial conditions of the crater shape, i.e. 
$R^*_0(1)$ (figure 15d) and 
$R^*_1(1)$ (figure 15e). Besides, the coefficient 
$R^*_0$ is the least modified with a maximum change of 
$\sim 30\,\%$ (figure 15iv), while 
$\phi ^*_0$, 
$\phi ^*_1$, 
$\phi ^*_2$ and 
$R^*_1$ reach respectively 
${\sim }300\,\%$ (figure 15i), 
${\sim }150\,\%$ (figure 15ii), 
$\sim 200\,\%$ (figure 15iii) and 
${\sim }100\,\%$ (figure 15v). Finally, the change in the coefficients over time is not homogeneous. For example, 
$\phi ^*_0$ is changed relatively uniformly over time (in magnitude), independently of the modified initial condition, while 
$\phi ^*_2$ is changed much more heterogeneously and depends on the modified initial condition.
Open access
