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:

(2.1a–c)\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$:

(2.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:

(2.3)\begin{align} u_r(r, \theta, t)&=\sum_{l=0}^{+\infty} u_{r,l}(r, t) \bar{P}_l(\cos\theta), \end{align}

(2.4)\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.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):

(3.1)\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

(3.2)\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.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

(4.1)\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

(4.2)\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

(4.3)\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

(4.4)\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:

(4.5)\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

(4.6)\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:

(4.7)\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.