2 Experimental methods

In our experiments the impacting droplet is composed of either water or ethanol mixed with food dye (mass fraction ${<}2\,\%$ ) for visualization purposes. The radius of the droplet $R_{d}$ is fixed to 0.9 mm for ethanol. The radius of the water droplets is fixed to 1.4 mm for most experiments and to 1.75 mm occasionally. The impacting droplet is released from a nozzle above the substrate. Small oscillations may be induced during the pinch-off of the droplet. Nevertheless the effect of these oscillations may safely be neglected in the inertia-driven spreading regime we focus on in our experiments. The impact speed, $U$ , reaches from $1.5~\text{m}~\text{s}^{-1}$ to $5.5~\text{m}~\text{s}^{-1}$ by altering the falling height. The granular target consists of ceramic beads with the specific density $\unicode[STIX]{x1D70C}_{g}=6000~\text{kg}~\text{m}^{-3}$ . Three different grain sizes are used, and the hydrophilicity of the grains is enhanced after cleaned with a piranha solution (see table 1) (Zhao, de Jong & van der Meer Reference Zhao, de Jong and van der Meer2017). The packing density of the bed $\unicode[STIX]{x1D719}$ is prepared in the range of 0.59–0.63 by air fluidization and subsequent mechanical taps. While the droplet deformation is visualized with a high-speed camera, at the same instance, the deformation of the target surface is measured by an in-house built high-speed laser profilometer (Zhao et al. Reference Zhao, de Jong and van der Meer2015b ). The impacting droplet appears to be a dark ellipse (due to perspective) in the video. When detecting the spreading diameter cross-correlations are computed between dark ellipse templates and the (background subtracted) image, and the template size corresponding to the maximum correlation gives the spreading diameter. The evolution of the crater profile with time, $t$ , is quantified as a function $z(r,t)$ . It is good to stress that here, as well as in all of the analyses in this study, axisymmetry of the impact crater is assumed. Here, $z$ represents the vertical coordinate, and $r$ is the in-plane radial coordinate with the impact centre at $r=0$ . When the location of impact is far from the laser lines, the part of $z(r,t)$ close to the impact centre is unresolved. We discard the experiments where the unresolved region is beyond $r=1~\text{mm}$ , and extrapolate the profile to $r=0$ with a parabolic fit otherwise, where a parabola centred at $r=0$ is fitted to the last 20 datapoints closest to the centre. The parabolic fitting function is given in appendix C. This fit is implemented on the available data in the range of $r<2~\text{mm}$ . To estimate the potential error the same method is applied to a fully resolved experimental profile, and the difference between the fitted crater depth and the measured one is very small (0.0043 mm), and is representative of potential errors made by using the fit. Note that the corresponding error in the crater area is even smaller, since the fitted region of $r<1~\text{mm}$ , the crater centre, is typically smaller than the region of interest (see below).

Table 1. Summary of the contact angle of water and ethanol, $\unicode[STIX]{x1D703}_{w}$ and $\unicode[STIX]{x1D703}_{e}$ , and the grain size, $d_{g}$ .

Figure 2. The leftmost column (a–c) illustrates the crater profile at the moment of maximal droplet expansion $t=t^{\ast }$ (dashed lines) and the moment the final static state of the crater is reached, $t=t^{\,f}$ (solid lines). For (a,b), the differences of crater profiles at these two moments result from the deformation of the liquid–grain mixtures, whereas the discrepancy at the crater rim in (c) is caused by grain splashing. The middle column (d–f) shows the time evolution of crater area, $\unicode[STIX]{x1D6F4}$ (see text for definition). The moment $t^{\ast }$ of maximum droplet expansion is indicated by vertical dashed lines in (d–f). In (e) the dash dotted line denotes $\unicode[STIX]{x1D6F4}(t^{\,f})$ . Both columns correspond to the impacts in figure 1. The rightmost column (g–i) displays vertical cross-sections from X-ray tomograms through the impact centre where the residual is highlighted with yellow colour. These tomograms are scanned for droplet impacts on glass beads, different experiments from that in (a–f). Details can be found in appendix A.

Figure 1 exemplifies three residual morphologies: truffle, pie and pancake. Those residual shapes result from different post-impact processes, as can be seen in the experimental movies in the supplementary materials available online at https://doi.org/10.1017/jfm.2019.742. With the acquired crater profile, $z(r,t)$ , differences between those processes can be quantified. We first compare the profiles at two moments in time. One moment is at $t=t^{\ast }$ , where the impacting droplet reaches its maximum expansion radius $R_{d}^{\ast }$ . The other moment is the last frame of the experiment, $t=t^{\,f}$ , where the crater has obtained its final static state. Minutes or hours afterwards the experiment, the liquid–grain residual is dried. Its shape may only change slightly from that at $t^{\,f}$ . The crater profiles of experiments shown in figure 1 are presented in figure 2(a–c). While the profiles at $t^{\ast }$ and $t^{\,f}$ are largely identical for the pancake shape (cf. figure 2 c), both truffle and pie shapes display prominent evolution (cf. figure 2 a,b). The difference between the truffle and pie shapes is clear as well. For the truffle shape the evolution occurs through the whole crater, from rim to the centre, and the final profile is highly curved. In contrast, for the pie shape, the part of the profile far from the impact centre remains pinned, and the final profile is relatively flat, with the centre lifted upwards.

To further quantify the crater evolution we calculate the reduced crater area by integrating the crater profile between $r=0$ and $r=R_{0}$ , $\unicode[STIX]{x1D6F4}(t)=\unicode[STIX]{x03C0}\int _{0}^{R_{0}}\sqrt{1+z^{\prime 2}}\,\text{d}r^{2}-\unicode[STIX]{x03C0}R_{0}^{2}$ . Here, $z^{\prime }=\unicode[STIX]{x2202}z(r,t)/\unicode[STIX]{x2202}r$ is the derivative of the profile with respect to $r$ , and $R_{0}$ is the radial distance where the crater profile intersects with $z=0$ at $t=t^{\ast }$ (cf. figure 2 b). The results corresponding to figure 2(a–c) are given in (d–f) in the same figure. As the flat area, $A_{0}=\unicode[STIX]{x03C0}R_{0}^{2}$ , is subtracted, $\unicode[STIX]{x1D6F4}(t)$ always increases from 0 upon impact to a local maximum at $t=t^{\ast }$ for all three cases (cf. figure 2 d–f). Afterwards the time evolution $\unicode[STIX]{x1D6F4}(t)$ exhibits interestingly distinct signatures corresponding to the three typical crater profile evolutions described above. The profile of pancake shape does not change between $t^{\ast }$ and $t^{\,f}$ , therefore, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}=\unicode[STIX]{x1D6F4}(t^{\,f})-\unicode[STIX]{x1D6F4}(t^{\ast })=0$ (cf. figure 2 f). In contrast, for the truffle shape, the liquid–grain mixture contracts into a highly curved residual, more curved than the crater profile at $t^{\ast }$ , which results in a positive area difference, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}>0$ (cf. figure 2 d). For the pie shape, as the centre of the crater is lifted, the profile becomes flatter, and $\unicode[STIX]{x1D6F4}$ decreases. Therefore, the area difference is negative, i.e. $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}<0$ (cf. figure 2 e). The lift of the crater centre creates a cavity revealed by the X-ray tomography scan after $t^{\,f}$ in figure 2(h) (refer to appendix A for details of the X-ray tomogram). Hence, the sign of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ can be used to indicate morphology: $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}=0$ for the pancake shape, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}>0$ for the truffle shape and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}<0$ for the pie shape. In experiments the oscillation of the remaining liquid surface introduces fluctuations into $\unicode[STIX]{x1D6F4}(t)$ after $t>t^{\ast }$ . When computing $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ we average the data of $\unicode[STIX]{x1D6F4}(t)$ in the last 2 ms to obtain the value $\unicode[STIX]{x1D6F4}(t^{\,f})$ at $t^{\,f}$ .

Figures 1 and 2 indicate that the final morphology is the outcome of the post-impact dynamics of the liquid–grain mixture. What leads to the different post-impact dynamics? From the formation of truffle morphology it is clear that in this case the dynamics is dominated by the Laplace pressure, the product of the surface tension and local curvature of the liquid–air interface, which are developed during the impact. This also give us clues to the formation of the other two morphologies.

3 Mixing ratio

Similar to impacts on a solid surface, upon impact the droplet starts to spread, and the circumference curvature is developed along its edges. Meanwhile, a central curvature also arises, which is unique to impacts on deformable substrates, like the granular target studied here (cf. figure 3 inset). These two curvatures determine the post-impact dynamics. According to the corresponding crater dimensions (Zhao et al. Reference Zhao, de Jong and van der Meer2015b , Reference Zhao, de Jong and van der Meer2017; de Jong, Zhao & van der Meer Reference de Jong, Zhao and van der Meer2017), the circumference curvature is always larger than the central one. Nevertheless, their relative significance can be altered by the degree of mixing between liquid and grains during the impact. When there is little mixing, the circumference curvature is dominant and leads to retraction of the liquid–grain mixture, which finally results in a ball-like shape, a truffle. In the other extreme of complete mixing, where the whole droplet penetrates into the granular substrate during the impact, both curvatures are lost, and no further dynamics occurs after the impact. The final morphology preserves the maximum expanded shape, the so-called pancake shape. However, in the intermediate range of mixing, where more than half of the droplet is mixed with grains, the circumference curvature is suppressed, which prevents the retraction. The radial extension of the residual is then frozen in the expanded state. Nevertheless, the central curvature of the top liquid–air interface may still lift the liquid–grain mixture and create the underneath cavity displayed in figure 2(h), the pie morphology. In this dual-curvature model the morphology appears to be a function of mixing ratio at $t^{\ast }$ , namely the ratio of the liquid volume that is mixed with grains, ${\mathcal{V}}_{m}$ , to that of the whole droplet, $V_{d}=4\unicode[STIX]{x03C0}R_{d}^{3}/3$ .

Figure 3. Main plot: crater area difference, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ , plotted versus the estimated mixing ratio ${\mathcal{V}}_{m}/V_{d}$ . The sign of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ , here normalized by the corresponding flat area $A_{0}$ , indicates the final morphology. The morphology develops with the estimated mixing ratio ${\mathcal{V}}_{m}/V_{d}$ . Three regimes can be distinguished which are in accordance with the residue morphology shown in figure 1: a, truffle; b, pie; c, pancake. See text for details and for the definition of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ and $A_{0}$ . Inset: a sketch of the two curvatures occurring in a spread droplet. In practice, the droplet (blue) is mixed with grains which are not drawn here.

To estimate ${\mathcal{V}}_{m}$ , we first decompose it into the contact area, ${\mathcal{A}}_{c}$ , and the penetration depth of the liquid into the substrate, $L$ . Although in reality the contact area evolves with time in a complicated manner, for the purpose of this analysis it is approximated by a constant, namely the geometric mean of the initial droplet radius $R_{d}$ and its maximum spreading radius $R_{d}^{\ast }$ , i.e. ${\mathcal{A}}_{c}=\unicode[STIX]{x03C0}R_{d}R_{d}^{\ast }$ . The development of $L$ is modelled by Darcy’s law (refer to appendix B), and its solution is

(3.1) $$\begin{eqnarray}L(t)=\sqrt{\frac{2\unicode[STIX]{x1D705}P}{\unicode[STIX]{x1D707}_{l}(1-\unicode[STIX]{x1D719})}\,t}.\end{eqnarray}$$

In the above equation the permeability of the substrate, $\unicode[STIX]{x1D705}=(1-\unicode[STIX]{x1D719})^{3}d_{g}^{2}/(180\unicode[STIX]{x1D719}^{2})$ , is defined by the Carman–Kozeny relation (Carman Reference Carman1956), and $\unicode[STIX]{x1D707}_{l}$ is the dynamic viscosity of the liquid. The driving pressure $P$ consists of two parts: the inertia pressure and capillary pressure, where the effect of wettability of grains on the liquid penetration is included in the latter. The complete expression of $P$ can be found in appendix B and Zhao et al. (Reference Zhao, de Jong and van der Meer2017). In the model the spreading time is estimated as half of the intrinsic capillary oscillation time of the droplet (Richard, Clanet & Quere Reference Richard, Clanet and Quere2002; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003; Delon et al. Reference Delon, Terwagne, Dorbolo, Vandewalle and Caps2011), $\unicode[STIX]{x1D70F}_{c}={\textstyle \frac{1}{2}}\sqrt{(\unicode[STIX]{x03C0}/6)(\unicode[STIX]{x1D70C}_{l}D_{0}^{3}/\unicode[STIX]{x1D70E})}$ , where $\unicode[STIX]{x1D70E}$ is the surface tension of the liquid. With all the quantities in (3.1) defined, we could finally evaluate ${\mathcal{V}}_{m}=(1-\unicode[STIX]{x1D719}){\mathcal{A}}_{c}L(\unicode[STIX]{x1D70F}_{c})$ . Note that the factor $(1-\unicode[STIX]{x1D719})$ accounts for the presence of grains in the mixture. The resultant ${\mathcal{V}}_{m}$ increases with impact speed $U$ , grain size $d_{g}$ and wettability, $\cos \unicode[STIX]{x1D703}_{c}$ .

The normalized crater area development $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}/A_{0}$ is plotted against mixing ratio ${\mathcal{V}}_{m}/V_{d}$ in figure 3. Three regimes corresponding to individual morphologies can be defined according to the sign of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}$ : a. $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}>0$ , the truffle regime; b. $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}<0$ , the pie regime; and c. $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}\approx 0$ , the pancake regime. As morphology develops with the mixing ratio, the first cross-over, between regime a and b, is found at ${\mathcal{V}}_{m}/V_{d}\approx 0.55$ , where approximately half of the droplet is mixed with grains, and the circumference curvature is lost. As a result, surface tension is not able to retract the droplet anymore, as discussed in Zhao et al. (Reference Zhao, de Jong and van der Meer2015b ). The second cross-over, between b and c, is at ${\mathcal{V}}_{m}/V_{d}\approx 1$ , where the whole droplet penetrates into the granular substrate, and the liquid–grain mixture is completely frozen at its maximum expansion. Because of mass conservation, ${\mathcal{V}}_{m}/V_{d}>1$ is physically clearly not expected, and it simply indicates that the whole droplet penetrates into the granular target before reaching $t=\unicode[STIX]{x1D70F}_{c}$ . In fact, this is consistent with a transition of the droplet spreading dynamics. It has been shown that the viscous dissipation of liquid penetrating into the granular target becomes the largest energy sink in this regime (Zhao et al. Reference Zhao, de Jong and van der Meer2017). Hence, regime c may also be interpreted from an energy perspective: no further dynamics would happen after the impact energy is largely dissipated by the viscosity of the liquid moving in the pores between the grains.

4 Lifting criterion

Both the wettability and the size of the grains play a role in the degree of mixing, and thus in the formation of final morphologies. In general, increasing the grain size or the wettability enhances the mixing ratio. For instance, impacts on $d_{g}=257~\unicode[STIX]{x03BC}\text{m}$ result in the pancake shape, while impacts on $d_{g}=98~\unicode[STIX]{x03BC}\text{m}$ do not enter that regime. Although the dual-curvature model introduced above captures this general tendency of the morphology development, in regime b ( $0.55<{\mathcal{V}}_{m}/V_{0}<1$ ) there are impacts of relatively large $d_{g}$ and low $U$ that actually end in the pancake shape ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}\approx 0$ ) rather than the expected pie morphology. Therefore, the mixing ratio alone is not sufficient to explain the formation of pie morphology. The missing factor is the force balance. If the crater is too shallow, the Laplace pressure offered by the central curvature cannot overcome the gravitational pressure, then no lifting would occur.

The lifting of the liquid–grain mixture is a process where surface energy, $E_{s}$ is converted into gravitational energy $E_{g}$ . The condition under which lifting can be initiated is given by $\text{d}(E_{s}+E_{g})/\text{d}t|_{t=t^{\ast }}<0$ . This condition is path dependent. To obtain a criterion in terms of accessible experimental quantities we introduce a virtual lifting process, where the concave liquid–grain mixture at $t^{\ast }$ is raised to be flat. Following this deformation the centre depth of the crater $Z_{c}$ is decreased from $Z_{c}^{\ast }$ to 0. Note that here depth is used for the downwards direction. Using the chain rule and $\text{d}Z_{c}/\text{d}t<0$ the lifting criterion can be rewritten now as

(4.1) $$\begin{eqnarray}\left.\frac{\text{d}(E_{s}+E_{g})}{\text{d}Z_{c}}\right|_{Z_{c}=Z_{c}^{\ast }}>0.\end{eqnarray}$$

If the crater profile is parabolic, $z(r,Z_{c})=Z_{c}(r^{2}/R_{0}^{2}-1)$ , one can readily obtain the proportionalities $E_{g}\propto V(Z_{c})\propto Z_{c}$ and $E_{s}\propto \unicode[STIX]{x1D6F4}(Z_{c})\propto Z_{c}^{2}$ (details in appendix C). Furthermore, the excess surface energy $E_{s}$ decreases from $\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6F4}^{\ast }$ to 0, and the gravitational energy $E_{g}$ increases from 0 to $\unicode[STIX]{x1D70C}_{m}\bar{L}gV^{\ast }$ (cf. figure 4). Here, $g=9.8~\text{m}~\text{s}^{-2}$ is the gravitational acceleration, $\unicode[STIX]{x1D70C}_{m}=\unicode[STIX]{x1D70C}_{l}(1-\unicode[STIX]{x1D719})+\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x1D719}$ is estimated as the density of the mixture and $\bar{L}={\mathcal{V}}_{m}/[\unicode[STIX]{x03C0}{R_{d}^{\ast }}^{2}(1-\unicode[STIX]{x1D719})]=LR_{d}/R_{d}^{\ast }$ is the average mixing layer thickness. Also, $\unicode[STIX]{x1D6F4}^{\ast }=\unicode[STIX]{x1D6F4}(t^{\ast })$ and $V^{\ast }=\unicode[STIX]{x03C0}\int _{0}^{R_{0}}z(r,t^{\ast })\,\text{d}r^{2}$ are the crater area and the crater volume at $t^{\ast }$ or when $Z_{c}=Z_{c}^{\ast }$ . Let $\unicode[STIX]{x0394}E_{s}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6F4}^{\ast }$ and $\unicode[STIX]{x0394}E_{g}=\unicode[STIX]{x1D70C}_{m}\bar{L}gV^{\ast }$ denote the total virtual change of $E_{s}$ and $E_{g}$ respectively. Both energies can then be written explicitly as a function of $Z_{c}$ $E_{g}=\unicode[STIX]{x0394}E_{g}(1-Z_{c}/Z_{c}^{\ast })$ and $E_{s}=\unicode[STIX]{x0394}E_{s}(Z_{c}/Z_{c}^{\ast })^{2}$ . Substituting this into (4.1) results in the lifting criterion,

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x0394}E_{s}}{\unicode[STIX]{x0394}E_{g}}>0.5.\end{eqnarray}$$

Note that $\unicode[STIX]{x0394}E_{s}$ and $\unicode[STIX]{x0394}E_{g}$ can be computed from the measured crater profile at $t^{\ast }$ .

Figure 4. A sketch of the virtual lifting process. The lifting process is denoted by the decrease of $Z_{c}$ . The lifting may occur if the total energy $E_{s}+E_{g}$ decreases at the beginning of the lifting, i.e. equation (4.1).

Figure 5. The data in regimes b and c in figure 3 are plotted versus the ratio of the virtual change of the surface energy, $\unicode[STIX]{x0394}E_{s}$ , and the gravitational energy, $\unicode[STIX]{x0394}E_{g}$ . See text for their definitions. The same symbols and colours as that in figure 3 are used here. The dashed line indicates the predicted lifting criterion, $\unicode[STIX]{x0394}E_{s}/\unicode[STIX]{x0394}E_{g}=0.5$ , which separates the pie and pancake shapes in the overlapping range of mixing ratio.

The measured energy ratio of the data in regime b and regime c in figure 3 is plotted against the mixing ratio in figure 5. It can be seen that the data points corresponding to the pancake shape in regime b ( $0.55<{\mathcal{V}}_{m}/V_{d}<1$ ) are at or below the proposed lifting criterion $\unicode[STIX]{x0394}E_{s}/\unicode[STIX]{x0394}E_{g}=0.5$ , whereas the pie data all lie well above. For a given mixing ratio the energy ratio is proportional to the product of two impact parameters: the area-to-volume ratio, $\unicode[STIX]{x1D6F4}^{\ast }/V^{\ast }$ , and squared the maximum expansion radius of the droplet, ${R_{d}^{\ast }}^{2}$ , i.e. $\unicode[STIX]{x0394}E_{s}/\unicode[STIX]{x0394}E_{g}\propto \unicode[STIX]{x1D6F4}^{\ast }{R_{d}^{\ast }}^{2}/V^{\ast }$ . Both $\unicode[STIX]{x1D6F4}^{\ast }/V^{\ast }$ and $R_{d}^{\ast }$ increase with impact speed $U$ (refer to appendix C). Therefore, the impacts with large $d_{g}$ and/or higher wettability but small $U$ in regime b fail to meet the lifting criterion.

With the established criterion we could discuss its possible application to the raindrop impression mentioned in the beginning. In nature, raindrops fall at their terminal velocity. For a droplet radius of 1.4 mm, the terminal velocity is approximately $7~\text{m}~\text{s}^{-1}$ . First, we evaluate the mixing ratio as a function of grain size $d_{g}$ using (3.1). Parameters used in the evaluation are the packing density of the target $\unicode[STIX]{x1D719}=0.59$ , the contact angle of the grains $\cos \unicode[STIX]{x1D703}_{c}=0.4$ and the maximum spreading radius of the droplet $R_{d}^{\ast }/R_{d}=4$ . It is possible to reach regime b of the mixing ratio for the grain size range of $35{-}75~\unicode[STIX]{x03BC}\text{m}$ . Silt falls in such a size range (Boggs Reference Boggs1995), and its mineral origin is quartz and feldspar, lighter than ceramic beads used here. Therefore, the density of the mixture $\unicode[STIX]{x1D70C}_{m}$ would be smaller. Furthermore, the larger impact speed results in a larger $R_{d}^{\ast }$ , and the area volume ratio, $\unicode[STIX]{x1D6F4}^{\ast }/V^{\ast }$ would be comparable to the impacts shown here (appendix C). It is therefore very likely that the resultant $\unicode[STIX]{x0394}E_{s}/\unicode[STIX]{x0394}E_{g}$ would have satisfied the lift criterion. In conclusion, it is possible for raindrop impact on a substrate like silt to result in the pie morphology with a underneath cavity. It is however worth pointing out that the curvature of the final pie shape is unlikely to exceed that of the liquid–grain mixture at $t^{\ast }$ (although the sign is opposite). The radius of curvature of the final morphology is larger than both the horizontal and vertical dimensions, $R_{d}^{\ast }$ and $Z_{c}^{\ast }$ . Therefore, the underneath cavity created by droplet impacts should be well distinguishable from that resulting from a nearly spherical bubble in the mud.