1. Introduction
Drop impacts have piqued the interest of scientists and artists alike for centuries, with the phenomenon being sketched by da Vinci (Reference da Vinci1508) in the early 16th century and photographed by Worthington (Reference Worthington1876a,Reference Worthingtonb) in the late 19th century. It is, indeed, captivating to observe raindrops hitting a solid surface (Kim et al. Reference Kim, Wu, Esmaili, Dombroskie and Jung2020; Lohse & Villermaux Reference Lohse and Villermaux2020) or ocean spray affecting maritime structures (Berny et al. Reference Berny, Popinet, Séon and Deike2021; Villermaux, Wang & Deike Reference Villermaux, Wang and Deike2022). The phenomenology of drop impact is extremely rich, encompassing behaviours such as drop deformation (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Moláček & Bush Reference Moláček and Bush2012), spreading (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Wildeman et al. Reference Wildeman, Visser, Sun and Lohse2016), splashing (Xu, Zhang & Nagel Reference Xu, Zhang and Nagel2005; Riboux & Gordillo Reference Riboux and Gordillo2014; Thoraval et al. Reference Thoraval, Schubert, Karpitschka, Chanana, Boyer, Sandoval-Naval, Dijksman, Snoeijer and Lohse2021), fragmentation (Villermaux & Bossa Reference Villermaux and Bossa2011; Villermaux Reference Villermaux2020), bouncing (Richard & Quéré Reference Richard and Quéré2000; Kolinski, Mahadevan & Rubinstein Reference Kolinski, Mahadevan and Rubinstein2014; Chubynsky et al. Reference Chubynsky, Belousov, Lockerby and Sprittles2020; Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020; Sharma & Dixit Reference Sharma and Dixit2021; Sanjay, Chantelot & Lohse Reference Sanjay, Chantelot and Lohse2023a) and wetting (de Gennes Reference de Gennes1985; Fukai et al. Reference Fukai, Shiiba, Yamamoto, Miyatake, Poulikakos, Megaridis and Zhao1995; Quéré Reference Quéré2008; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). These behaviours are influenced by the interplay of inertial, capillary and viscous forces, as well as additional factors like non-Newtonian properties (Bartolo, Josserand & Bonn Reference Bartolo, Josserand and Bonn2005; Bartolo et al. Reference Bartolo, Boudaoud, Narcy and Bonn2007; Smith & Bertola Reference Smith and Bertola2010; Gorin et al. Reference Gorin, Di Mauro, Bonn and Kellay2022) of the liquid and even ambient air pressure (Xu et al. Reference Xu, Zhang and Nagel2005), making the parameter space for this phenomenon both extensive and high-dimensional.
Naturally, even the process of a Newtonian liquid drop impacting a rigid substrate is governed by a plethora of control parameters, including but not limited to the drop's density 
$\rho _d$, diameter 
$D_0$, velocity 
$V_0$, dynamic viscosity 
$\eta _d$, surface tension 
$\gamma$ and acceleration due to gravity 
$g$ (figure 1a). To navigate this rich landscape, we focus on two main dimensionless numbers that serve as control parameters (figure 1b): the Weber number 
$We$, which is the ratio of inertial to capillary forces and is given by
\begin{equation} We = \frac{\rho_d V_0^2 D_0}{\gamma}; \end{equation}
and the Ohnesorge number 
$Oh$, which captures the interplay between viscous damping and capillary oscillations, offering insights into how viscosity affects the drop's behaviour upon impact,
\begin{equation} Oh = \frac{\eta_d}{\sqrt{\rho_d \gamma D_0}}. \end{equation}Additionally, the Bond number
\begin{equation} Bo = \frac{\rho_dgD_0^2}{\gamma} \end{equation}compares gravity with inertial forces and is needed to uniquely define the non-dimensional problem.
(a) Problem schematic with an axisymmetric computational domain used to study the impact of a drop with diameter 
$D_0$ and velocity 
$V_0$ on a non-wetting substrate. In the experiments, we use a quartz force sensor to measure the temporal variation of the impact force. The subscripts 
$d$ and 
$a$ denote the drop and air, respectively, to distinguish their material properties, which are the density 
$\rho$ and the dynamic viscosity 
$\eta$. The drop–air surface tension coefficient is 
$\gamma$. The grey dashed–dotted line represents the axis of symmetry, 
$r = 0$. Boundary air outflow is applied at the top and side boundaries (tangential stresses, normal velocity gradient and ambient pressure are set to zero). The domain boundaries are far enough from the drop not to influence its impact process (
$\mathcal {L}_{max} \gg D_0$, 
$\mathcal {L}_{max} = 8R$ in the worst case). (b) The phase space with control parameters: the Weber number (
$We$, dimensionless kinetic energy) and the Ohnesorge number (
$Oh$, dimensionless viscosity), exemplifying different applications.
The drop impact is not only interesting from the point of view of fundamental research but also finds relevance in inkjet printing (Lohse Reference Lohse2022), the spread of respiratory drops carrying airborne microbes (Bourouiba Reference Bourouiba2021; Ji, Yang & Feng Reference Ji, Yang and Feng2021; Pöhlker et al. Reference Pöhlker2023), cooling applications (Kim Reference Kim2007; Shiri & Bird Reference Shiri and Bird2017; Jowkar & Morad Reference Jowkar and Morad2019), agriculture (Bergeron et al. Reference Bergeron, Bonn, Martin and Vovelle2000; Bartolo et al. Reference Bartolo, Boudaoud, Narcy and Bonn2007; Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018; Sijs & Bonn Reference Sijs and Bonn2020; He et al. Reference He, Ding, Li, Mu, Li and Liu2021; Hoffman et al. Reference Hoffman, Sijs, de Goede and Bonn2021), criminal forensics (Smith & Brutin Reference Smith and Brutin2018; Smith, Nicloux & Brutin Reference Smith, Nicloux and Brutin2018) and many other industrial and natural processes (Rein Reference Rein1993; Yarin Reference Yarin2006; Tuteja et al. Reference Tuteja, Choi, Ma, Mabry, Mazzella, Rutledge, McKinley and Cohen2007; Cho et al. Reference Cho, Preston, Zhu and Wang2016; Hao et al. Reference Hao, Liu, Chen, Li, Zhang, Zhao and Wang2016; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Liu, Wang & Jiang Reference Liu, Wang and Jiang2017; Yarin, Roisman & Tropea 2017; Yarin et al. Reference Yarin, Roisman and Tropea2017; Wu et al. Reference Wu, Du, Alsaid, Wu, Hua, Yan, Yao, Ma, Zhu and He2020). For these applications, it is pertinent to understand the forces involved in drop impacts, as these forces can lead to soil erosion (Nearing, Bradford & Holtz Reference Nearing, Bradford and Holtz1986) or damage to engineered surfaces (Gohardani Reference Gohardani2011; Ahmad, Schatz & Casey Reference Ahmad, Schatz and Casey2013; Amirzadeh et al. Reference Amirzadeh, Louhghalam, Raessi and Tootkaboni2017). We refer the readers to Cheng, Sun & Gordillo (Reference Cheng, Sun and Gordillo2022) for an overview of the recent studies unravelling drop impact forces (see also Li et al. (Reference Li, Zhang, Guo and Lv2014), Soto et al. (Reference Soto, de Lariviere, Boutillon, Clanet and Quéré2014), Philippi, Lagrée & Antkowiak (Reference Philippi, Lagrée and Antkowiak2016), Zhang et al. (Reference Zhang, Li, Guo and Lv2017), Gordillo, Sun & Cheng (Reference Gordillo, Sun and Cheng2018), Mitchell et al. (Reference Mitchell, Klewicki, Korkolis and Kinsey2019) and Zhang et al. (Reference Zhang, Zhang, Lv, Li and Guo2019)).
These forces have been studied by Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022), employing experiments and simulations and deriving scaling laws. A liquid drop impacting a non-wetting substrate undergoes a series of phases – spreading, recoiling and potentially rebounding (Chantelot Reference Chantelot2018) – driven by the normal reaction force exerted by the substrate (figure 2). The moment of touchdown (see figure 2a,b, 
$t = 0$ to 
$t_1$) (Wagner Reference Wagner1932; Philippi et al. Reference Philippi, Lagrée and Antkowiak2016; Gordillo et al. Reference Gordillo, Sun and Cheng2018) is not surprisingly associated with a pronounced peak in the temporal evolution of the drop impact force 
$F(t)$ owing to the sudden deceleration as high as 
$100$ times the acceleration due to gravity (Clanet et al. Reference Clanet, Béguin, Richard and Quéré2004) (figure 2a, 
$F_1/(\rho _dV_0^2D_0^2) \approx 0.82, 0.92, 0.99$ for 
$Oh = 0.0025, 0.06\ \text {and}\ 0.2$, respectively; at 
$t_1 \approx 0.03\sqrt {\rho _dD_0^3/\gamma }$). The force diminishes as the drop reaches its maximum spreading diameter (figure 2a,b, 
$t = t_m$). Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) revealed that also the jump-off is accompanied by a peak in the normal reaction force, which was up to then unknown (figure 2a, 
$F_2/(\rho _dV_0^2D_0^2) \approx 0.37, 0.337, 0.1$ for 
$Oh = 0.0025, 0.06\ \text {and}\ 0.2$, respectively; for the second force peak amplitude – at time 
$t_2 \approx 0.42\sqrt {\rho _dD_0^3/\gamma }$ after impact). The second peak in the force also coincides with the formation of a Worthington jet, a narrow upward jet of liquid that can form due to flow focusing by the retracting drop (figure 2a,b, 
$t = t_2$). Under certain conditions (
$We \approx 9$, 
$Oh = 0.0025$), this peak can be even more pronounced than the first. This discovery is critical for superhydrophobicity which is volatile and can fail due to external disturbances such as pressure (Lafuma & Quéré Reference Lafuma and Quéré2003; Callies & Quéré Reference Callies and Quéré2005; Sbragaglia et al. Reference Sbragaglia, Peters, Pirat, Borkent, Lammertink, Wessling and Lohse2007; Li et al. Reference Li, Quéré, Lv and Zheng2017), evaporation (Tsai et al. Reference Tsai, Lammertink, Wessling and Lohse2010; Chen et al. Reference Chen, Ma, Li, Hao, Guo, Luk, Li, Yao and Wang2012; Papadopoulos et al. Reference Papadopoulos, Mammen, Deng, Vollmer and Butt2013), mechanical vibration (Bormashenko et al. Reference Bormashenko, Pogreb, Whyman and Erlich2007) or the impact forces of prior droplets (Bartolo et al. Reference Bartolo, Bouamrirene, Verneuil, Buguin, Silberzan and Moulinet2006a).
Comparison of the drop impact force 
$F(t)$ obtained from experiments and simulations for the three typical cases with impact velocity 
$V_0 = 1.2, 0.97, 0.96\ {\rm m}\ {\rm s}^{-1}$, diameter 
$D_0 = 2.05, 2.52, 2.54\ {\rm mm}$, surface tension 
$\gamma \approx 72, 61, 61\ {\rm mN}\ {\rm m}^{-1}$ and viscosity 
$\eta _d = 1, 25.3, 80.2\ {\rm mPa}\ {\rm s}$. These parameters give 
$Oh = 0.0025, 0.06, 0.2$ and 
$We = 40$. For the three cases, the two peak amplitudes, 
$F_1/(\rho _dV_0^2D_0^2) \approx 0.82, 0.92, 0.99$ at 
$t_1 \approx 0.03\sqrt {\rho _dD_0^3/\gamma }$ and 
$F_2/(\rho _dV_0^2D_0^2) \approx 0.37, 0.337, 0.1$ at 
$t_2 \approx 0.42\sqrt {\rho _dD_0^3/\gamma }$, characterize the inertial shock from impact and the Worthington jet before takeoff, respectively. The drop reaches the maximum spreading at 
$t_{max}$ when it momentarily stops and retracts until 
$0.8\sqrt {\rho _dD_0^3/\gamma }$ when the drop takes off (
$F = 0$). The black and grey dashed lines in (a) mark 
$F = 0$ and the resolution 
$F = 0.5\ {\rm mN}$ of our piezoelectric force transducer, respectively. (b) Four instances are further elaborated through numerical simulations for (
$We = 40, Oh = 0.0025$), namely (i) 
$t = 0$ (touchdown), (ii) 
$t = 0.03\sqrt {\rho _dD_0^3/\gamma }$ (
$t_1$), (iii) 
$t = 0.2\sqrt {\rho _dD_0^3/\gamma }$ (
$t_{max}$) and (iv) 
$t = 0.42\sqrt {\rho _dD_0^3/\gamma }$ (
$t_2$). The insets of (a) exemplify these four instances for the three representative cases illustrated here. The experimental snapshots are overlaid with the drop boundaries from simulations. We stress the excellent agreement between experiments and simulations without any free parameters. The left-hand part of each numerical snapshot shows (on a 
$\log _{10}$ scale) the dimensionless local viscous dissipation function 
$\tilde {\xi }_\eta \equiv \xi _\eta D_0/(\rho _dV_0^3) = 2Oh(\boldsymbol {\tilde {\mathcal {D}}:\tilde {\mathcal {D}}})$, where 
$\boldsymbol {\mathcal {D}}$ is the symmetric part of the velocity gradient tensor, and the right-hand part the velocity field magnitude normalized with the impact velocity. The black velocity vectors are plotted in the centre of mass reference frame of the drop to clearly elucidate the internal flow. Also see supplementary movies SM1–SM3 available at https://doi.org/10.1017/jfm.2024.982.
In contrast to our prior study Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022), which fixed the Ohnesorge number to that of a 2 mm diameter water drop (
$Oh = 0.0025$), our present investigation reported in this paper explores a broader parameter space. We systematically and independently vary the Weber and Ohnesorge numbers, extending the range of 
$Oh$ to as high as 100. This comprehensive approach enables us to develop new scaling laws and provides a more unified understanding of the forces involved in drop impact problems. Our findings are particularly relevant for applications with varying viscosities and impact velocities (figure 1).
The structure of this paper is as follows. Section 2 briefly describes the experimental and numerical methods. Sections 3 and 4 offer detailed analyses of the first and second peaks, respectively, focusing on their relationships with the Weber number (
$We$) and the Ohnesorge number (
$Oh$). Conclusions and perspectives for future research are presented in § 5.
2. Methods
2.1. Experimental method
In the experimental set-up, shown schematically in figure 1(a), a liquid drop impacts a superhydrophobic substrate. For water drops, such a surface is coated with silanized silica nanobeads with a diameter of 20 nm (Glaco Mirror Coat Zero; Soft99) resulting in the advancing and receding contact angles of 
$167 \pm 2^{\circ }$ and 
$154 \pm 2^{\circ }$, respectively (Gauthier et al. Reference Gauthier, Symon, Clanet and Quéré2015; Li et al. Reference Li, Quéré, Lv and Zheng2017). On the other hand, for viscous aqueous glycerine drops, the upper surface is coated with an acetone solution of hydrophobic beads (Ultra ever Dry, Ultratech International, a typical bead size of 20 nm), resulting in the advancing and receding contact angles of 
$166 \pm 4^{\circ }$ and 
$159 \pm 2^{\circ }$, respectively (Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). The properties of the impacting drop are controlled using water–glycerine mixtures with viscosities 
$\eta _d$ varying by almost two orders of magnitude, from 1 to 
$80.2\ {\rm mPa}\ {\rm s}$. Surface tension is either 
$72\ {\rm mN}\ {\rm m}^{-1}$ (pure water) or 
$61\ {\rm mN}\ {\rm m}^{-1}$ (glycerol), while density 
$\rho _d$ ranges from 1000 to 
$1220\ {\rm kg}\ {\rm m}^{-3}$, as detailed in table 1 (Cheng Reference Cheng2008; Volk & Kähler Reference Volk and Kähler2018; Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). We note that using liquids such as silicone oil can provide a broader range of viscosity variation when paired with a superamphiphobic substrate (Deng et al. Reference Deng, Mammen, Butt and Vollmer2012). Additionally, employing drops of smaller radii facilitates the exploration of higher Ohnesorge numbers (
$Oh$, see (1.2)). The drop diameter 
$D_0$ is controlled between 2.05 and 
$2.76\ {\rm mm}$ by pushing it through a calibrated needle (see Appendix A for details). Consequently, we calculate 
$Oh$ using the properties in table 1. The Weber number (
$We$, see (1.1)) is set using the impact velocity 
$V_0$ varying between 0.38 and 
$2.96\ {\rm m}\ {\rm s}^{-1}$ by changing the release height of the drops above the substrate. All experiments are conducted at ambient pressure and temperature. The impact force is directly measured using a high-precision piezoelectric force transducer (Kistler 9215A) with a resolution of 
$0.5\ {\rm mN}$. During these measurements, the high-frequency vibrations induced by the measurement system and the surrounding noise are spectrally removed using a low-pass filter with a cutoff frequency of 
$5\ {\rm kHz}$, following the procedure in Li et al. (Reference Li, Zhang, Guo and Lv2014), Zhang et al. (Reference Zhang, Li, Guo and Lv2017), Gordillo et al. (Reference Gordillo, Sun and Cheng2018) and Mitchell et al. (Reference Mitchell, Klewicki, Korkolis and Kinsey2019). The experiment also employs a high-speed camera (Photron Fastcam Nova S12) synchronized at 10 000 f.p.s. with a shutter speed 1/20 000 s. Throughout the manuscript, the error bars are of a statistical nature (one standard deviation) and originate from repeated trials. They are visible if they are larger than the marker size. We refer the readers to the supplementary material of Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) and Appendix A for further details of the experimental set-up and error characterization of the dimensionless control parameters.
Properties of the water–glycerol mixtures used in the experiments. Here 
$\rho _d$ and 
$\eta _d$ are the density and viscosity of the drop, respectively, and 
$\gamma$ denotes the liquid–air surface tension coefficient (Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). These properties are calculated using the protocol provided in Cheng (Reference Cheng2008) and Volk & Kähler (Reference Volk and Kähler2018).
2.2. Numerical framework
In the direct numerical simulations employed for this study, the continuity and the momentum equations take the form
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} = 0 \end{equation}and
\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{vv}) = \frac{1}{\rho}(-\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot}(2\eta{\boldsymbol{\mathcal{D}}}) + \boldsymbol{f}_\gamma) + \boldsymbol{g}, \end{equation}
respectively. Here, 
$\boldsymbol {v}$ is the velocity field, 
$t$ is time, 
$p$ is pressure and 
$\boldsymbol {g}$ is acceleration due to gravity. We use the free software program Basilisk C that employs the well-balanced geometric volume of fluid (VoF) method (Popinet Reference Popinet2009, Reference Popinet2018). The VoF tracer 
$\varPsi$ delineates the interface between the drop (subscript d, 
$\psi = 1$) and air (subscript a, 
$\psi = 0$), introducing a singular force 
$\boldsymbol {f}_\gamma \approx \gamma \kappa \boldsymbol {\nabla }\varPsi$ (
$\kappa$ denotes interfacial curvature, see Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992)) to respect the dynamic boundary condition at the interface. This VoF tracer sets the material properties such that density 
$\rho$ and viscosity 
$\eta$ are given by
\begin{equation} \rho = \rho_a + (\rho_d - \rho_a)\varPsi \end{equation}and
\begin{equation} \eta = \eta_a + (\eta_d - \eta_a)\varPsi, \end{equation}respectively. This VoF field is advected with the flow, following the equation
\begin{equation} \frac{\partial \varPsi}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{v}\varPsi) = 0. \end{equation}
Lastly, we calculate the normal reaction force 
$\boldsymbol {F}(t)$ by integrating the pressure field 
$p$ at the substrate,
\begin{equation} \boldsymbol{F}(t) = \left(\int_\mathcal{A} (\,p-p_0)\,{\rm d}\mathcal{A}\right)\hat{\boldsymbol{z}}, \end{equation}
where, 
$p_0$, 
$d\mathcal {A}$ and 
$\hat {\boldsymbol {z}}$ are the ambient pressure, substrate area element and the unit vector normal to the substrate, respectively.
We leverage the axial symmetry of the drop impact (figure 1a). This axial symmetry breaks at large 
$We$ (
${\ge }100$ for water drops and even larger Weber number for more viscous drops), owing to destabilization by the surrounding gas after splashing (Xu et al. Reference Xu, Zhang and Nagel2005; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Driscoll & Nagel Reference Driscoll and Nagel2011; Riboux & Gordillo Reference Riboux and Gordillo2014; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). To solve the governing equations ((2.1)–(2.5)), the velocity field 
$\boldsymbol {v}$ and time 
$t$ are normalized by the inertiocapillary scales, 
$V_\gamma = \sqrt{\gamma/(\rho_dD_0)}$ and 
$\tau _\gamma = \sqrt {\rho _dD_0^3/\gamma }$, respectively. Furthermore, the pressure is normalized using the capillary pressure scale 
$p_\gamma = \gamma /D_0$. In such a conceptualization, 
$Oh$ and 
$We$ described in § 1 uniquely determine the system. The Ohnesorge number based on air viscosity 
$Oh_a = (\eta _a/\eta _d)Oh$ and air–drop density ratio 
$\rho _a/\rho _d$ are fixed at 
$10^{-5}$ and 
$10^{-3}$, respectively, to minimize the influence of the surrounding medium on the impact forces. Lastly, we keep the Bond number 
$Bo$ (see (1.3)) fixed at 1 throughout the manuscript. In our system, the relevance of gravity is characterized by the dimensionless Froude number 
$Fr = V_0^2/(gD_0) = We/Bo$ which compares inertia with gravity. Throughout this manuscript, 
$Fr > 1$ and gravity's role is subdominant compared with inertia (for detailed discussion, see Appendix B). The substrate is modelled as a no-slip and non-penetrable wall, whereas vanishing stress and pressure are applied at the remaining boundaries to mimic outflow conditions for the surrounding air. The domain boundaries are far enough from the drop not to influence its impact process (
$\mathcal {L}_{max} \gg D_0$, 
$\mathcal {L}_{max} = 8R$ in the worst case). At 
$t = 0$, in our simulations, we release a spherical drop whose south pole is 
$0.05D_0$ away from the substrate and is falling with a velocity 
$V_0$. It is important to note that large experimental drops may deviate from perfect sphericity due to air drag as they fall after detaching from the needle and potential residual oscillations from detachment. These shape perturbations are more pronounced in cases with low Weber and Ohnesorge numbers. To quantify this non-sphericity, we measure the drop's aspect ratio (horizontal to vertical diameter) immediately before substrate contact. The precise preimpact drop shape can significantly influence subsequent impact dynamics (Thoraval et al. Reference Thoraval, Takehara, Etoh and Thoroddsen2013; Yun Reference Yun2017; Zhang et al. Reference Zhang, Zhang, Lv, Li and Guo2019). In our experiments, we constrain our analysis to drops with aspect ratios between 0.96 and 1.05. Given this narrow range, we posit that the impact of these shape variations is negligible compared with the experimental error bars derived from repeated trials under identical nominal conditions. The simulations utilize adaptive mesh refinement to finely resolve the velocity, viscous dissipation and the VoF tracer fields. A minimum grid size 
$\varDelta = D_0/2048$ is used for this study.
To ensure a perfectly non-wetting surface, we impose a thin air layer (minimum thickness 
${\sim }\varDelta /2$) between the drop and the substrate. This air layer prevents direct contact between the liquid and solid (Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014; Sprittles Reference Sprittles2024), effectively mimicking a perfectly non-wetting surface. The presence of this air layer is crucial for capturing the dynamics of drop impact on superhydrophobic surfaces, as it allows for the formation of an air cushion that can significantly affect the spreading and rebound behaviour of the drop (Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020; Sanjay et al. Reference Sanjay, Chantelot and Lohse2023a). While this approach does not fully resolve the microscopic dynamics within the air layer itself, such as the high-velocity gradients and viscous dissipation inside the gas film once it thins below a critical size (
${\sim }10\varDelta$), it has been shown to accurately capture the macroscopic behaviour of drop impact in the parameter range of interest (Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020; Alventosa, Cimpeanu & Harris Reference Alventosa, Cimpeanu and Harris2023; Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b; García-Geijo, Riboux & Gordillo Reference García-Geijo, Riboux and Gordillo2024). We refer the readers to Sanjay (Reference Sanjay2022) for discussions about this ‘precursor’ air film method and to Popinet Collaborators (Reference Popinet2013–2023), Sanjay (Reference Sanjay2024) and Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) for details on the numerical framework.
3. Anatomy of the first impact force peak
This section elucidates the anatomy of the first impact force peak and its relationship with the Weber 
$We$ and Ohnesorge 
$Oh$ numbers, first for the inertial limit (§ 3.1, 
$Oh \ll 1$) and then for the viscous asymptote (§ 3.2, 
$Oh \gg 1$). The results of this section are summarized in figure 3 that shows an excellent agreement between experiments and simulations without any free parameters.
Anatomy of the first impact force peak amplitude at low 
$Oh$ in between 0.0025 and 0.2, see colour legend: 
$We$ dependence of the (a) magnitude 
$F_1$ normalized by the inertial force scale 
$F_\rho = \rho _dV_0^2D_0^2$ and time 
$t_1$ to reach the first force peak amplitude normalized by (b) the inertial time scale 
$\tau _\rho = D_0/V_0$ and (c) the inertiocapillary time scale 
$\tau _\gamma = \sqrt {\rho _dD_0^3/\gamma }$. The Jupyter notebook for producing the figure can be found: https://www.cambridge.org/S0022112024009820/JFM-Notebooks/files/figure3/figure3.ipynb.
3.1. Low Ohnesorge number impacts
For low 
$Oh$ and large 
$We$, inertial force and time scales dictate the drop impact dynamics (figures 3 and 4). As the drop falls on a substrate, the part of the drop immediately in contact with the substrate stops moving, whereas the top of the drop still falls with the impact velocity (figure 4, from 
$t = t_1/4$ until 
$t = t_1$). Consequently, momentum conservation implies
\begin{equation} F_1 \sim V_0\frac{{\rm d} m}{{\rm d} t}, \end{equation}
where the mass flux 
${\rm d} m/{\rm d} t \sim \rho _d V_0D_0^2$ (Soto et al. Reference Soto, de Lariviere, Boutillon, Clanet and Quéré2014; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). As a result, the first peak amplitude scales with the inertial pressure force (figure 3a)
\begin{equation} F_1 \sim F_\rho,\quad \text{where }F_\rho = \rho_d V_0^2D_0^2, \end{equation}
for high Weber numbers (
$We > 30$, 
$F_1 \approx 0.81F_\rho$). Furthermore, the time 
$t_1$ to reach 
$F_1$ follows
\begin{equation} t_1 \sim \frac{D_0}{V_0} = \tau_\rho, \end{equation}
where, 
$\tau _\rho$ is the inertial time scale. The relation between (3.2) and (3.3) is apparent from the momentum conservation which implies that the impulse of the first force peak is equal to the momentum of the impacting drop, i.e. 
$F_1t_1 \sim \rho _dV_0D_0^3 = F_\rho \tau _\rho$ (see Gordillo et al. (Reference Gordillo, Sun and Cheng2018), Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) and figure 3b,c). These scaling laws depend only on the inertial shock at impact and are wettability-independent (Zhang et al. Reference Zhang, Li, Guo and Lv2017; Gordillo et al. Reference Gordillo, Sun and Cheng2018; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). For details of the scaling law, including the prefactors, we refer the readers to Philippi et al. (Reference Philippi, Lagrée and Antkowiak2016), Gordillo et al. (Reference Gordillo, Sun and Cheng2018) and Cheng et al. (Reference Cheng, Sun and Gordillo2022).
Direct numerical simulations snapshots illustrating the drop impact dynamics for 
$We = (\textit {a})$ 40 and (b) 2. The left-hand side of each numerical snapshot shows the pressure normalized by (a) the inertial pressure scale 
$\rho _dV_0^2$ and (b) the capillary pressure scale 
$\gamma /D_0$. The right-hand side shows the velocity field magnitude normalized by the impact velocity 
$V_0$.
Figure 3 further illustrates that this inertial asymptote is insensitive to viscosity variations up to 100-fold as 
$F_1 \sim F_\rho$ and 
$t_1 \sim \tau _\rho$ for 
$0.0025 < Oh < 0.2$. However, deviations from the inertial force and time scales are apparent for 
$We < 30$ (figure 3), a phenomenon also reported in earlier work (Soto et al. Reference Soto, de Lariviere, Boutillon, Clanet and Quéré2014; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). In these instances, inertia does not act as the sole governing force but instead complements surface tension, which dictates the pressure inside the drop (
$\,p \sim \gamma /D_0$ throughout the drop for 
$We \lesssim 1$, figure 4b). Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) proposed an empirical functional dependence as
\begin{equation} F_1 = \left(\alpha_1\rho_d V_0^2 + \alpha_2\frac{\gamma}{D_0}\right)D_0^2, \end{equation}
based on dimensional analysis, with 
$\alpha _1$ and 
$\alpha _2$ as free parameters which were determined to be approximately 
$1.6$ and 
$0.81$, respectively, for water (
$Oh = 0.0025$). These coefficients only deviate marginally in the current work despite the significant increase in 
$Oh$ as compared with previous works (Cheng et al. Reference Cheng, Sun and Gordillo2022; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). This consistency underscores the invariance of the pressure field inside the drop to an increase in 
$Oh$ (close to the impact region, figure 4a and throughout the drop, figure 4b).
3.2. Large Ohnesorge number impacts
Figure 5 reaffirms the findings of § 3.1 for low 
$Oh$ that the first impact peak amplitude 
$F_1$ and the time to reach this peak amplitude 
$t_1$ scale with 
$F_\rho$ and 
$\tau _\rho$, respectively. As the Ohnesorge number increases further, the first impact force peak amplitude normalized with 
$F_\rho$ begins to increase, indicating a transition around 
$Oh \approx 0.1$, where viscosity starts to play a significant role. At large 
$Oh$, we observe the scaling relationship (figure 5a)
\begin{equation} F_1 \sim F_\rho\sqrt{Oh}. \end{equation}
The drop's momentum is still 
$\rho _dV_0D_0^3$ which must be balanced by the impulse from the substrate, 
$F_1t_1$ (see § 3.1, Gordillo et al. (Reference Gordillo, Sun and Cheng2018), Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) and figure 5b). Consequently, the time 
$t_1$ follows
\begin{equation} t_1 \sim \frac{\tau_\rho}{\sqrt{Oh}}. \end{equation}Anatomy of the first impact force peak amplitude for viscous impacts from our numerical simulations: the 
$Oh$ dependence of (a) the magnitude 
$F_1$ normalized by the inertial force scale 
$\rho _dV_0^2D_0^2$ and (b) the time 
$t_1$ to reach the first force peak amplitude normalized by inertial time scale 
$\tau _\rho = D_0/V_0$. (c) The 
$Re$ dependence of the magnitude 
$F_1$ normalized by the inertial force scale 
$\rho _dV_0^2D_0^2$ as compared with the (implicit) theoretical calculation of Gordillo et al. (Reference Gordillo, Sun and Cheng2018). The black line corresponds to the scaling relationship described in § 3.2. The Weber number is colour-coded.
Figure 5 further shows that these scaling laws are weakly dependent on the Weber number, as viscous dissipation consumes the entire initial kinetic energy of the impacting drop (figure 6). Once again, we stress that using the water–glycerol mixtures limits the range of 
$Oh$ that we can probe experimentally. We further note that the first peak is robust and does not depend on the wettability of the substrate. Consequently, to compare with the existing data such as those in Cheng et al. (Reference Cheng, Sun and Gordillo2022) with different liquids to cover a wider range of liquid viscosities and to account for the apparent 
$We$-dependence, we plot 
$F_1$ compensated with 
$F_\rho$ against the impact Reynolds number 
$Re \equiv \sqrt {We}/Oh = V_0D_0/\nu _d$. For the low 
$Re$ regime, such a plot allows us to describe the 
$We$ dependence on the prefactor more effectively, as illustrated in figure 5(c). However, it is important to note that some scatter is still observed at high 
$Re$ values, which can be attributed to the 
$We$ dependence of the impact force peak amplitude. This lack of a pure scaling behaviour demonstrates how the interplay between kinetic energy and viscous dissipation within the drop dictates the functional dependence of the maximum impact force on 
$Oh$.
Direct numerical simulations snapshots illustrating the drop impact dynamics for 
$We = 100$ and 
$Oh = (\textit {a})$ 0.05, (b) 0.5 and (c) 5. The left-hand side of each numerical snapshot shows the viscous dissipation function 
$\xi _\eta$ normalized by the inertial scale 
$\rho _dV_0^3/D_0$. The right-hand side shows the velocity field magnitude normalized by the impact velocity 
$V_0$.
To systematically elucidate these scaling behaviours in the limit of small 
$Re$, we need to find the typical scales for the rate of change of kinetic energy and that of the rate of viscous dissipation for the drop impact system. First, we can readily define an average rate of viscous dissipation per unit mass as
\begin{equation} \bar{\varepsilon} \sim \frac{1}{\tau_\rho}\frac{1}{D_0^3} \int_0^{\tau_\rho}\int_\varOmega\nu_d(\boldsymbol{\mathcal{D}:\mathcal{D}})\,{\rm d}\varOmega \,{\rm d} t, \end{equation}
where 
$\nu _d$ is the kinematic viscosity of the drop and 
$d\varOmega$ is the volume element where dissipation occurs. Notice that 
$\bar {\varepsilon }$ has the dimensions of 
$V_0^3/D_0$, i.e. length squared over time cubed or velocity squared over time, as it should be for dissipation rate of energy per unit mass. We can estimate 
$\varOmega = D_{foot}^2l_\nu$ (figure 6), where 
$D_{foot}$ is the drop's foot diameter in contact with the substrate and 
$l_\nu$ is the viscous boundary layer thickness. This boundary layer marks the region of strong velocity gradients (
${\sim }V_0/l_\nu$) analogous to the Mirels (Reference Mirels1955) shockwave-induced boundary layer. For details, we refer the authors to Schlichting (Reference Schlichting1968), Schroll et al. (Reference Schroll, Josserand, Zaleski and Zhang2010) and Philippi et al. (Reference Philippi, Lagrée and Antkowiak2016). Consequently, the viscous dissipation rate scales as
\begin{equation} \bar{\varepsilon} \sim \frac{1}{\tau_\rho D_0^3}\int_0^{\tau_\rho}\nu_d \left(\frac{V_0}{l_\nu}\right)^2 D_{foot}^2l_\nu \,{\rm d} t. \end{equation}
To calculate 
$D_{foot}$, we assume that the drop maintains a spherical cap shape throughout the impact (figure 6). To calculate the distance the drop would have travelled if there were no substrate, we use the relation 
$d \sim V_0t$. Simple geometric arguments allow us to determine the relation between the foot diameter and this distance, 
$D_{foot} \sim \sqrt {D_0d}$ (Lesser Reference Lesser1981; Mandre, Mani & Brenner Reference Mandre, Mani and Brenner2009; Zheng, Dillavou & Kolinski Reference Zheng, Dillavou and Kolinski2021; Bertin Reference Bertin2024; Bilotto et al. Reference Bilotto, Kolinski, Lecampion, Molinari, Subhash and Garcia-Suarez2023). Interestingly, this scaling behaviour is similar to the inertial limit (Wagner Reference Wagner1932; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Philippi et al. Reference Philippi, Lagrée and Antkowiak2016; Gordillo, Riboux & Quintero Reference Gordillo, Riboux and Quintero2019) as discussed by Langley, Li & Thoroddsen (Reference Langley, Li and Thoroddsen2017) and Bilotto et al. (Reference Bilotto, Kolinski, Lecampion, Molinari, Subhash and Garcia-Suarez2023). Furthermore, the viscous boundary layer 
$l_\nu$ can be approximated using 
$\sqrt {\nu _d t}$ (Mirels Reference Mirels1955; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Philippi et al. Reference Philippi, Lagrée and Antkowiak2016). Filling these in (3.8), we get
\begin{equation} \bar{\varepsilon} \sim \frac{1}{\tau_\rho D_0^2}\int_0^{\tau_\rho}\sqrt{\nu_d} V_0^3 \sqrt{t} \,{\rm d} t, \end{equation}which on integration gives
\begin{equation} \bar{\varepsilon} \sim \sqrt{\nu_d \tau_\rho}V_0^3/D_0^2, \end{equation}
where 
$\tau _\rho$ is the inertial time scale. Here, we assume that for highly viscous drops, all energy is dissipated within a fraction of 
$\tau _\rho$. Filling in (3.10) and normalizing 
$\bar {\varepsilon }$ with the inertial scales 
$V_0^3/D_0$,
\begin{equation} \frac{\bar{\varepsilon}}{V_0^3/D_0} \sim \sqrt{\frac{\nu_d\tau_\rho}{D_0^2}} = \frac{1}{\sqrt{Re}} = \left(\frac{Oh}{\sqrt{We}}\right)^{1/2}. \end{equation}Next, the kinetic energy of the falling drop is given by
\begin{equation} \dot{K}(t) \equiv \frac{{\rm d} K(t)}{{\rm d} t} \sim \rho_dD_0^3\bar{\varepsilon},\quad\text{where } K(t) = \frac{1}{2}m(V(t))^2, \end{equation}
and 
$V(t)$ is the drop's centre of mass velocity. The left-hand side of (3.12) can be written as
\begin{equation} \dot{K}(t) = mV(t)\frac{{\rm d} V(t)}{{\rm d} t} = F(t)V(t). \end{equation}
In (3.13), 
$F(t)$ and 
$V(t)$ scale with the first impact force peak amplitude 
$F_1$ and the impact velocity 
$V_0$, respectively, giving the typical scale of the rate of change of kinetic energy as
\begin{equation} \dot{K}^* \sim F_1V_0. \end{equation}
We stress that (3.14) states that the rate of change of kinetic energy is equal to the power of the normal reaction force, an observation already made by Wagner (Reference Wagner1932) and Philippi et al. (Reference Philippi, Lagrée and Antkowiak2016) in the context of impact problems. Lastly, at large 
$Oh$, viscous dissipation enervates kinetic energy completely giving (figure 6c, also see Philippi et al. (Reference Philippi, Lagrée and Antkowiak2016) and Wildeman et al. (Reference Wildeman, Visser, Sun and Lohse2016))
\begin{equation} \dot{K}^* \sim F_1V_0 \sim \rho_dD_0^3\bar{\varepsilon}. \end{equation}Additionally, we use the inertial scales to non-dimensionalize (3.15) and fill in (3.11), giving
\begin{equation} \frac{F}{F_\rho} \sim \frac{\bar{\varepsilon}}{V_0^3/D_0} \sim \frac{1}{\sqrt{Re}} = \left(\frac{Oh}{\sqrt{We}}\right)^{1/2} \end{equation}
and using 
$F_1t_1 \sim \rho _dV_0D_0^3 = F_\rho \tau _\rho$,
\begin{equation} \frac{t_1}{\tau_\rho} \sim \left(\frac{\sqrt{We}}{Oh}\right)^{1/2}. \end{equation}In summary, we use energy and momentum invariance to elucidate the parameter dependencies of the impact force as illustrated in figure 5. The scaling arguments capture the dominant force balance during the impact process, considering the relative importance of inertial, capillary and viscous forces. As the dimensionless viscosity of impacting drops increases, the lack of surface deformation increases the normal reaction force (3.16). Further, the invariance of incoming drop momentum implies that this increase in normal reaction force occurs on a shorter time scale (3.17).
4. Anatomy of the second impact force peak
This section delves into the anatomy of the second impact force peak amplitude 
$F_2$ as a function of the Weber 
$We$ and Ohnesorge 
$Oh$ numbers, summarized in figure 7. We once again note the remarkable agreement between experiments and numerical simulations in this figure.
Anatomy of the second impact force peak amplitude: 
$We$ dependence of the (a) magnitude 
$F_2$ normalized by the inertial force scale 
$F_\rho = \rho _dV_0^2D_0^2$ and time 
$t_2$ to reach the second force peak amplitude normalized by (b) the inertiocapillary time scale 
$\tau _\gamma = \sqrt {\rho _dD_0^3/\gamma }$ and (c) inertial time scale 
$\tau _\rho = D_0/V_0$. The Jupyter notebook for producing the figure can be found: https://www.cambridge.org/S0022112024009820/JFM-Notebooks/files/figure7/figure7.ipynb.
Similar to the mechanism leading to the formation of the first peak (§ 3), also the mechanism for the formation of this second peak is momentum conservation. As the drop takes off from the surface, it applies a force on the substrate. As noted in § 1 and Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022), this force also coincides with the formation of a Worthington jet (figure 2 iv). The time 
$t_2$ at which the second peak is observed scales with the inertiocapillary time scale and is insensitive to 
$We$ and 
$Oh$ (figure 7b,c). Once again, we invoke the analogy between drop oscillation and drop impact to explain this behaviour (Richard, Clanet & Quéré Reference Richard, Clanet and Quéré2002; Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012). At the time instant 
$t_2 \approx 0.44\tau _\gamma$, the drop's internal motion undergoes a transition from a predominantly radial flow to a vertical one due to the formation of the Worthington jet (Chantelot Reference Chantelot2018; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022). Figure 8 exemplifies this jet in the 
$Oh$–
$We$ parameter space, which is intricately related to the second peak in the drop impact force. For low 
$Oh$ and large 
$We$, the drop retraction follows a modified Taylor–Culick dynamics (Bartolo et al. Reference Bartolo, Josserand and Bonn2005; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Sanjay et al. Reference Sanjay, Sen, Kant and Lohse2022). As 
$We$ is increased, the jet gets thinner but faster, maintaining a constant momentum flux 
$\rho _dV_j^2d_j^2$, where 
$V_j$ and 
$d_j$ are the jet's velocity and diameter, respectively (figure 8, Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022)). This invariance leads to the observed scaling 
$F_2 \sim F_\rho$ in this regime (
$F_2 \approx 0.37F_\rho$ for 
$We \ge 30, Oh \le 0.01$).
Direct numerical simulations snapshots illustrating the influence of 
$We$ and 
$Oh$ on the inception of the Worthington jet. All these snapshots are taken at the instant when the second peak appears in the temporal evolution of the normal reaction force (
$t = t_2$). The left-hand side of each numerical snapshot shows the viscous dissipation function 
$\xi _\eta$ normalized by the inertial scale 
$\rho _dV_0^3/D_0$. The right-hand side shows the velocity field magnitude normalized by the impact velocity 
$V_0$. The grey velocity vectors are plotted in the centre of mass reference frame of the drop to clearly elucidate the internal flow.
Furthermore, the low 
$We$ and 
$Oh$ regime relies entirely on capillary pressure (figure 4). Subsequently, 
$F_2 \sim F_\gamma = \gamma D_0$ for 
$Oh < 0.01$ and 
$We < 30$ (figure 7). This flow focusing (figure 8) is most efficient for 
$We = 9$ (figure 9a, 
$t_2/2 < t < t_2$, Renardy et al. (Reference Renardy, Popinet, Duchemin, Renardy, Zaleski, Josserand, Drumright-Clarke, Richard, Clanet and Quéré2003) and Bartolo, Josserand & Bonn (Reference Bartolo, Josserand and Bonn2006b)) where the capillary resonance leads to a thin-fast jet, accompanied by a bubble entrainment, reminiscent of the hydrodynamic singularity (figure 9, Sanjay, Lohse & Jalaal (Reference Sanjay, Lohse and Jalaal2021) and Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022)). The characteristic feature of this converging flow is a higher magnitude of 
$F_2$ compared with 
$F_1$ (figure 7).
Direct numerical simulations snapshots illustrating the influence of 
$We$ and 
$Oh$ on the singular Worthington jet: (a) 
$(We, Oh) = (9, 0.0025)$; (b) 
$Oh = 0.0025$ with 
$We =$ (i) 
$5$ and (ii) 
$12$; (c) 
$We = 9$ with 
$Oh =$ (i) 
$0.005$ and (ii) 
$Oh = 0.05$. The left-hand side of each numerical snapshot shows the viscous dissipation function 
$\xi _\eta$ normalized by inertial scale 
$\rho _dV_0^3/D_0$. The right-hand side shows the velocity field magnitude normalized by the impact velocity 
$V_0$.
However, this singular jet regime is very narrow in the 
$Oh$–
$We$ phase space. Figure 9(b) shows two cases for water drops (
$Oh = 0.0025$) at different 
$We$ (
$We = 5$ and 
$12$ for figures 9b i and 9b ii, respectively). Bubble entrainment does not occur in either of these cases. Consequently, the maximum force amplitude diminishes for these two cases (figure 7). Nonetheless, these cases are still associated with high local viscous dissipation near the axis of symmetry owing to the singular nature of the flow. Another mechanism to inhibit this singular Worthington jet is viscous dissipation in the bulk. As the Ohnesorge number increases, this singular jet formation disappears (
$Oh = 0.005$, figure 9c i), significantly reducing the second peak of the impact force. For even higher viscosities, the drop no longer exhibits the sharp, focused jet formation seen at lower viscosities, and the second peak in the force is notably diminished (
$Oh = 0.05$, figure 9c ii).
Lastly, as 
$Oh$ increases, bulk dissipation becomes dominant (apparent from increasing 
$Oh$ at fixed 
$We$ in figure 8) and can entirely inhibit drop bouncing. Recently, Jha et al. (Reference Jha, Chantelot, Clanet and Quéré2020) and Sanjay et al. (Reference Sanjay, Chantelot and Lohse2023a) showed that there exists a critical 
$Oh$, two orders of magnitude higher than that of a 
$2\ {\rm mm}$ diameter water drop, beyond which drops do not bounce either, irrespective of their impact velocity. Consequently, the second peak in the impact force diminishes for larger 
$Oh$, which explains the monotonic decrease of the amplitude 
$F_2$ observed in figure 7 for 
$We > 30, Oh > 0.01$.
5. Conclusion and outlook
In this work, we study the forces and dissipation encountered during the drop impact process by employing experiments, numerical simulations and theoretical scaling laws. We vary the two dimensionless control parameters – the Weber (
$We$, dimensionless impact kinetic energy) and the Ohnesorge number (
$Oh$, dimensionless viscosity) independently to elucidate the intricate interplay between inertia, viscosity and surface tension in governing the forces exerted by a liquid drop upon impact on a non-wetting substrate.
For the first impact force peak amplitude 
$F_1$, owing to the momentum balance after the inertial shock at impact, figure 10(a) summarizes the different regimes in the 
$Oh$–
$We$ phase space. For low 
$Oh$, inertial forces predominantly dictate the impact dynamics, such that 
$F_1$ scales with the inertial force 
$F_\rho$ (Philippi et al. Reference Philippi, Lagrée and Antkowiak2016; Gordillo et al. Reference Gordillo, Sun and Cheng2018; Mitchell et al. Reference Mitchell, Klewicki, Korkolis and Kinsey2019; Cheng et al. Reference Cheng, Sun and Gordillo2022; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022) and is insensitive to viscosity variations up to 100-fold. As 
$Oh$ increases, the viscosity becomes significant, leading to a new scaling law: 
$F_1 \sim F_\rho \sqrt {Oh}$. The paper unravels this viscous scaling behaviour by accounting for the loss of initial kinetic energy owing to viscous dissipation inside the drop. Lastly, at low 
$We$, the capillary pressure inside the drop leads to the scaling 
$F_1 \sim F_\gamma$ (Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Moláček & Bush Reference Moláček and Bush2012).
Regime map in terms of the drop Ohnesorge number 
$Oh$ and the impact Weber number 
$We$ to summarize the two peaks in the impact force by showing the different regimes described in this work based on (a) the first peak in the impact force peak amplitude 
$F_1$ and (b) the second peak in the impact force peak amplitude 
$F_2$. Both peaks are normalized by the inertial force scale 
$F_\rho = \rho _dV_0^2D_0^2$. These regime maps are constructed using 
$\sim 1500$ simulations in the range 
$0.001 \leq Oh \leq 100$ and 
$1 \leq We \leq 1000$. The grey solid line in (a) and dashed line in (b) mark the inertial–viscous transition (
$Re = 1$) and the bouncing–no-bouncing transition (
$Oh_c = 0.53$ for 
$Bo = 1$, see Sanjay et al. (Reference Sanjay, Chantelot and Lohse2023a)), respectively.
The normal reaction force described in this work is responsible for deforming the drop as it spreads onto the substrate, where it stops thanks to surface tension. If the substrate is non-wetting, it retracts to minimize the surface energy and finally takes off (Richard & Quéré Reference Richard and Quéré2000). In this case, the momentum conservation leads to the formation of a Worthington jet and a second peak in the normal reaction force, as summarized in figure 10(b). For low 
$Oh$ and high 
$We$, the second force peak amplitude scales with the inertial force (
$F_\rho$), following a modified Taylor–Culick dynamics (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010). In contrast, capillary forces dominate at low 
$We$ and low 
$Oh$, leading to a force amplitude scaling of 
$F_2 \sim F_\gamma$. We also identify a narrow regime in the 
$Oh$–
$We$ phase space where a singular Worthington jet forms, significantly increasing 
$F_2$ (Bartolo et al. Reference Bartolo, Josserand and Bonn2006b; Zhang et al. Reference Zhang, Sanjay, Shi, Zhao, Lv and Lohse2022), localized in the parameter space for 
$We \approx 9$ and 
$Oh < 0.01$. As 
$Oh$ increases, bulk viscous dissipation counteracts this jet formation, diminishing the second peak and ultimately inhibiting drop bouncing.
Our findings have far-reaching implications, not only enriching the fundamental understanding of fluid dynamics of drop impact but also informing practical applications in diverse fields such as inkjet printing, public health, agriculture and material science where the entire range of 
$Oh$–
$We$ phase space is relevant (figures 1b and 10). While this has identified new scaling laws, it also opens avenues for future research. For instance, it would be interesting to use the energy accounting approach to unify the scaling laws for the maximum spreading diameter for arbitrary 
$Oh$ (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Wildeman et al. Reference Wildeman, Visser, Sun and Lohse2016). Although, the implicit theoretical model summarized in Cheng et al. (Reference Cheng, Sun and Gordillo2022) describes most of data in figure 5, we stress the importance of having a predictive model to determine 
$F_1$ for given 
$We$ and 
$Oh$ (Sanjay & Lohse Reference Sanjay and Lohse2024). The 
$We$ influence on the impact force also warrants further exploration, especially in the regime 
$We \ll 1$ for arbitrary 
$Oh$ (Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Moláček & Bush Reference Moláček and Bush2012) and drop impact on compliant surfaces (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023; Ma & Huang Reference Ma and Huang2023). Another potential extension of this work is to non-Newtonian fluids (Martouzet et al. Reference Martouzet, Jørgensen, Pelet, Biance and Barentin2021; Agüero et al. Reference Agüero, Alventosa, Harris and Galeano-Rios2022; Bertin Reference Bertin2024; Jin et al. Reference Jin, Zhao, Blin, Allais, Mouterde and Quéré2023).
Supplementary material
Computational Notebooks files are available as supplementary material at https://doi.org/10.1017/jfm.2024.982 and online at https://www.cambridge.org/S0022112024009820/JFM-Notebooks.
Acknowledgements
We would like to thank V. Bertin for comments on the manuscript. We would also like to thank A. Prosperetti, P. Chantelot, D. van der Meer and U. Sen for illuminating discussions.
Funding
We acknowledge the funding by the ERC Advanced Grant no. 740479-DDD and NWO-Canon grant FIP-II. B.Z. and C.L. are grateful for the support from National Natural Science Foundation of China (grant nos. 12172189, 11921002, 11902179). This work was carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research. This work was sponsored by NWO – Domain Science for the use of supercomputer facilities.
Code availability
The codes used in the present article, and the parameters and data to reproduce figures 3 and 7, are permanently available at Sanjay (Reference Sanjay2024).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Note on the error characterization for the control parameters
This appendix outlines the methodology for characterizing experimental errors in quantification of the drop's size and impact velocities which is crucial for accurate calculation of the dimensionless control parameters 
$We$ and 
$Oh$. The drop diameter determination involves multiple steps. First, we measure the total mass (
$M_{100}$) of 100 drops using an electric balance. From this mass, using the liquid density and assuming spherical shape, we calculated the drop diameter (
$D_0$). We repeated this process five times, yielding 
$D_{0,1}$ to 
$D_{0,5}$. The average of these measurements provided the final drop diameter (
$D_0$) and its standard error. For impact velocity determination, we extracted data from experimental high-speed imagery. By tracking the drop centre's position in successive frames prior to substrate contact, and knowing the frame rate, we calculated the impact velocity. We repeated this process for five trials, obtaining 
$V_{0,1}$ to 
$V_{0,5}$. The average of these values gave the final impact velocity (
$V_0$) and its standard error.
The standard errors for drop diameters do not exceed 
$0.13\ {\rm mm}$. For instance, drops with Ohnesorge numbers of 
$0.0025$, 
$0.06$ and 
$0.2$ have diameters of 
$2.05 \pm 0.13\ {\rm mm}$, 
$2.52 \pm 0.11\ {\rm mm}$ and 
$2.54 \pm 0.09\ {\rm mm}$, respectively. The standard errors for impact velocities did not exceed 
$0.02\ {\rm m}\ {\rm s}^{-1}$. For the same 
$Oh$ values, the impact velocities were 
$1.2 \pm 0.002\ {\rm m}\ {\rm s}^{-1}$, 
$0.97 \pm 0.01\ {\rm m}\ {\rm s}^{-1}$ and 
$0.96 \pm 0.01\ {\rm m}\ {\rm s}^{-1}$, respectively. The combined errors in 
$D_0$ and 
$V_0$ resulted in approximately 
${\pm }7\,\%$ error in Weber number 
$We$ and 
${\pm }3\,\%$ error in Ohnesorge number 
$Oh$. Consequently, the horizontal error bars, which relate to errors in the control parameters, are smaller than the symbol sizes in our figures.
Appendix B. Role of gravity on drop impact forces
Following table 1 and considering the variation in impacting drop diameter (Appendix A), the Bond number (1.3) in our experiments ranges from 
$0.5$ to 
$1.25$, introducing an additional dimensionless control parameter alongside 
$We$ and 
$Oh$. Gravity typically plays a negligible role in these impact processes (Sanjay et al. Reference Sanjay, Chantelot and Lohse2023a; Sanjay & Lohse Reference Sanjay and Lohse2024). We undertook a sensitivity test varying the Bond number from 
$0$ to 
$1.25$ in our simulations. Figure 11 confirms the leading-order Bond invariance of the results as the impact force profiles, including both force peaks 
$F_1$ and 
$F_2$ and their corresponding times 
$t_1$ and 
$t_2$, remain invariant to these Bond number variations. Notably, while gravity does play a role in drop impact dynamics, particularly for longer time scales and in determining the critical Ohnesorge number 
$Oh_c$ for bouncing inhibition (see figure 10b and Sanjay et al. (Reference Sanjay, Chantelot and Lohse2023a)), its effect on the initial impact force peaks is minimal for the parameter range studied here (large Froude numbers, 
$Fr > 1$). This Bond number invariance allows us to focus on the more dominant effects of Weber and Ohnesorge numbers. Consequently, we selected the representative value of 
$Bo = 1$, corresponding to a diameter of 
$0.00254\ {\rm mm}$, density 
$1000\ {\rm kg}\ {\rm m}^{-3}$, gravitational acceleration 
$10\ {\rm m}\ {\rm s}^{-2}$ and surface tension 
$0.061\ {\rm N}\ {\rm m}^{-1}$.
Comparison of the drop impact force 
$F(t)$ obtained from simulations for the four different Bond numbers 
$Bo = 0, 0.5, 1, 1.25$. Here, 
$Oh = 0.06$ and 
$We = 40$. Both force peaks 
$F_1$ and 
$F_2$ as well as time to reach these peaks 
$t_1$ and 
$t_2$ are invariant to variation in 
$Bo$.
Open access
