2.1. Linear quasi-potential model for the liquid bath

We assume viscous effects are small except in the lubrication layer, and follow methodology from Dias, Dyachenko & Zakharov (Reference Dias, Dyachenko and Zakharov2008) which we briefly outline below. Take a Helmholtz decomposition of the liquid as a small perturbation from potential flow $\boldsymbol {u}_b = \boldsymbol {\nabla } \phi _b + \boldsymbol {\nabla } \times \boldsymbol {\psi _b}$, and use tangential stress conditions and a boundary-layer argument to eliminate the $\boldsymbol {\psi }_b$ terms, expressing the rotational part of the flow in terms of $\phi _b$ and $\eta _b$. A more detailed explanation of this argument can be found in Phillips et al. (Reference Phillips, Cimpeanu and Milewski2024). The free surface of the bath is given by $z=\eta _b(\boldsymbol {x},t)$ where $z=0$ is its undisturbed free surface. The governing linear system then becomes

(2.5)$$\begin{gather} 0= \Delta_H \phi_b + \partial_z^2 \phi_b, \quad z\le0, \end{gather}$$

(2.6)$$\begin{gather}\partial_t \phi_b ={-}\frac{1}{\rho_b}p_a - g\eta_b + \frac{\sigma_b}{\rho_b}\Delta_H \eta_b + 2\nu_b \Delta_H\phi_b, \quad z=0, \end{gather}$$

(2.7)$$\begin{gather}\partial_t \eta_b = \partial_z \phi_b + 2\nu_b \Delta_H \eta_b, \quad z=0, \end{gather}$$

where $\Delta _H$ is used to denote the horizontal Laplacian, the curvature is small and of the form $\kappa _b = -\boldsymbol {\nabla }_H \eta _b$, and velocity continuity $\boldsymbol {u}_b=\boldsymbol {u}_a$ is imposed at the interface $\eta _b$ between the bath and air layer. The terms proportional to $\nu _b$ arise from $\boldsymbol {\psi _b}$ and are the corrections to the pressure and vertical velocity, respectively.

2.2. Air-layer lubrication approximation

Away from the droplet–bath interaction, the air is assumed to be atmospheric pressure and have negligible effect on the system, as we must introduce the domain in which the lubrication effects may be important. The lower part of the droplet is described by $z=\eta _d(\boldsymbol {x}_H,t) = \boldsymbol {X}(t)-S(\boldsymbol {x}_H,t)$, composed of the shape of the lower part of the droplet $S$, and the vertical height of its centre of mass $\boldsymbol {X}$. The domain of $S$ in $\boldsymbol {x}_H$ is chosen such that $S$ is single valued and with bounded gradients. The dynamics of the air layer are described through a lubrication approximation as in Phillips et al. (Reference Phillips, Cimpeanu and Milewski2024). Consider a lubrication region $\varOmega _a = \varOmega _L\times [\eta _b(r,\theta,t),\eta _d(r,\theta,t)]$ where its footprint is described in cylindrical coordinates as

(2.8)\begin{equation} \varOmega_{L} = \{(r,\theta) : r\in [0, r^*(\theta,t)], \theta \in [0,2{\rm \pi}) \}, \end{equation}

where $r^*$ is the edge of the lubrication region. Denoting $h=\eta _d-\eta _b$, the edge of the lubrication layer is given by a criterion $h(r^*,\theta,t) = \varepsilon$. It is assumed that $\varOmega _L$ is a subset of the domain of $S$. Assuming that the layer between the droplet and the bath is thin, the leading-order balance of the Navier–Stokes equations results in the lubrication balance where $\partial _z p_a = 0$ with

(2.9)$$\begin{gather} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{u}_a^H + \partial_z w_a = 0, \quad \boldsymbol{x}\in \varOmega_{a}, \end{gather}$$

(2.10)$$\begin{gather}\mu_a \partial_z^2 \boldsymbol{u}_a^H = \boldsymbol{\nabla}_H p_a, \quad \boldsymbol{x}\in \varOmega_{a}, \end{gather}$$

and velocity continuity holds at both interfaces: $\boldsymbol {u}_a = \boldsymbol {u}_d$ at $z =\eta _d$ and $\boldsymbol {u}_a = \boldsymbol {\nabla }_H \phi$ at $z =\eta _b$. Here $\boldsymbol {\nabla }_H$ is the horizontal gradient operator, and $\boldsymbol {u}_a^H = (u_r,u_\theta )$ is the horizontal velocity of the air layer, such that $\boldsymbol {u}_a = (\boldsymbol {u}_a^H,w_a)$. As expected, the pressure $p_a = P(r,\theta,t)$ is independent of $z$ throughout the air layer. Further, we can solve the equation for $\boldsymbol {u}_a^H$ in terms of $P$:

(2.11)\begin{equation} \boldsymbol{u}_a^H = \frac{\boldsymbol{\nabla}_H P}{2\mu_a}(z-\eta_b)(z-\eta_d) + \frac{\boldsymbol{\nabla}_H \phi_b|_{\eta_b}}{\eta_b- \eta_d}(z-\eta_d) + \frac{\boldsymbol{u}^H_d|_{\eta_d}}{\eta_d - \eta_b}(z-\eta_b). \end{equation}

To obtain the thin-film equation, we integrate vertically the conservation of mass (incompressibility) equation resulting in

(2.12) \begin{equation} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \int_{\eta_b} ^{\eta_d} \boldsymbol{u}_a^H \,\text{d}z +\Big(\boldsymbol{u}_a^H \Big|_{\eta_b} \boldsymbol{\cdot} \boldsymbol{\nabla}_H\eta_b -\boldsymbol{u}_a^H \Big|_{\eta_d} \boldsymbol{\cdot} \boldsymbol{\nabla}_H\eta_d \Big) + (w_a|_{\eta_d} - w_a|_{\eta_b})= 0. \end{equation}

All terms after the integral can be expressed exactly as $\partial _t\eta _d- \partial _t \eta _b = \partial _t h$ from the kinematic conditions at both interfaces. Finally, evaluating the integral in the first term using (2.11) to obtain the form of the flux

(2.13)\begin{equation} \boldsymbol{Q}= \int_{\eta_b}^{\eta_d} \boldsymbol{u}_a^H \,\text{d}z ={-}\frac{1}{12\mu_a}h^3 \boldsymbol{\nabla}_H P + \frac{1}{2}h (\boldsymbol{\nabla}_H\phi_b + \boldsymbol{u}_d^H) \end{equation}

results in the thin-film equation

(2.14)\begin{equation} \partial_t h + \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{Q} = 0, \quad \boldsymbol{x_H} \in \varOmega_{L}. \end{equation}

As we shall see below, this should be considered a nonlinear free-boundary elliptic equation for $P$ where $-h_t$ is considered as a forcing term, with the boundary condition $P=0$ at the evolving boundary curve $\partial \varOmega _{L}$.

2.3. Droplet model

2.3.1. Droplet deformation

Linear damped droplet oscillations are well understood (Lamb Reference Lamb1924) and the results have been used in a variety of contexts, e.g. Aalilija, Gandin & Hachem (Reference Aalilija, Gandin and Hachem2020). In this section we write equivalent equations for these oscillations using the same quasi-potential theory used above for the bath. Consider a capillary-scale droplet of liquid which when unperturbed is a sphere of radius $R_0$. Describing the problem using spherical polar coordinates with azimuthal angle $\theta$ and polar angle $\varphi$, we can take a Helmholtz decomposition for the liquid velocity $\boldsymbol {u}_d = \boldsymbol {\nabla } \phi _d + \boldsymbol {\nabla } \times \boldsymbol {\psi }_d$. Using the same argument argument as above (Phillips et al. Reference Phillips, Cimpeanu and Milewski2024) we obtain similar governing equations to those of the bath (2.5)–(2.7):

(2.15)$$\begin{gather} \Delta \phi_d = 0,\quad r< R_0, \end{gather}$$

(2.16)$$\begin{gather}\partial_t \phi_d ={-}\frac{1}{\rho_d}p_a - \frac{\sigma_d}{\rho_d} \kappa_d - 2\nu_d \partial_r^2 \phi_d,\quad r=R_0, \end{gather}$$

(2.17)$$\begin{gather}\partial_t \zeta_d = \partial_r \phi_d + 2\nu_d\left(\frac{1}{r^2}\Delta_s \zeta_d - \frac{1}{r^3}\int \Delta_s \phi_d \,\text{d}t\right), r=R_0, \end{gather}$$

where surface of the droplet is given by $r=R_0+\zeta _d(\theta,\varphi )$, and $\kappa _d$ is the mean curvature

(2.18a,b)\begin{equation} \kappa_d ={-} \left( \frac{2}{R_0^2} \zeta_d + \frac{1}{R_0^2} \Delta_s \zeta_d \right), \quad \Delta = \frac{1}{r^2}\partial_r (r^2 \partial_r) + \frac{1}{r^2}\Delta_s, \end{equation}

with $\Delta _s$ the Laplace–Beltrami operator (i.e. the surface Laplacian). The equations are valid in the small-oscillation limit and have been truncated to contain terms that are linear in the Ohnesorge number (see Appendix A for a discussion of the non-dimensional equations). One may proceed with an eigenfunction decomposition of $\zeta _d$ and $\phi _d$ as

(2.19a,b)\begin{equation} \phi_d = \sum_{l=2}^\infty \sum_{m={-}l}^l a_{lm}(t) r^l Y_l^m, \quad \zeta_d = \sum_{l=2}^\infty \sum_{m={-}l}^l c_{lm}(t) Y_l^m, \end{equation}

where the $a_{lm}(t), c_{lm}(t)$ are coefficients of the $Y_l^m (\theta,\varphi ) = {\rm e}^{{\rm i}m\varphi }P_l^m(\cos \theta )$ spherical harmonics that satisfy

(2.20)\begin{equation} {-}\Delta_s Y_l^m = (l)(l+1)Y_l^m. \end{equation}

The $P_l^m$ are the associated Legendre polynomials of degree $l$ and order $m$. The terms $l=0,1$ have been omitted as $l=0$ corresponds to a dilation of the sphere and $l=1$ corresponds to translational and rotational symmetries which need to be considered separately as they will be affected by the external forces due to the air layer. Substitution into the boundary conditions (2.16)–(2.17) leads to the differential equation (where leading order contributions have been determined with respect to the Ohnesorge number, Oh, and terms of order $Oh^2$ have been omitted)

(2.21)\begin{equation} \ddot{c}_{l,m} + 2\lambda_{l} \dot{c}_{l,m} + \omega_{l}^2c_{l,m} ={-}\frac{l\hat{p}_{l,m}}{\rho R_0}, \end{equation}

where $\hat {p}_{l,m}$ is the $Y_{l,m}$ coefficient of the lubrication air pressure $p_a$. The coefficients of the velocity potential are given to leading order by $a_{l,m}={\dot {c}_{l,m}}/{(lR_0^{l-1}})$. The damping coefficients and modal frequencies are given by

(2.22a,b)\begin{equation} \lambda_l = (2l+1)(l-1)\frac{\mu}{\rho R_0^2}, \quad \omega_l^2 = l(l-1)(l+2) \frac{\sigma}{\rho R_0^3}, \end{equation}

as seen in the literature (Lamb Reference Lamb1924; Aalilija et al. Reference Aalilija, Gandin and Hachem2020; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023). A similar equation may be obtained for the coefficients $a_{lm}(t)$; however, these are unnecessary for the lubrication-mediated (LM) model and have been omitted for brevity.

2.3.2. Equations for the centre of mass

The location of the droplet is determined through updating the position of the centre of mass, given by $\boldsymbol {X}(t)$. We shall disregard drag effects from the air, instead only considering contributions from the lubrication layer, then the equation of motion is

(2.23)\begin{equation} m\ddot{\boldsymbol{X}} = F ={-}\int_{\varOmega_L} p_a \, \boldsymbol{n_d} \,\text{d}a - mg \, \hat{\boldsymbol{z}}, \end{equation}

where $\boldsymbol {n_d}$ is the outward normal of the droplet, $m = 4\rho {\rm \pi}R_0^3/3$ is its mass, $g$ is the gravitational constant and $\hat {\boldsymbol {z}}$ is the unit vector in the vertical direction. The boundary of the droplet is now given by $\boldsymbol {X} + (R_0+\zeta _d(\theta,\varphi )) \boldsymbol {\hat {r}}$, where $\boldsymbol {\hat {r}}$ is a radial unit vector relative to the position $\boldsymbol {X}$.

We note that we have disregarded global rotational motion of the droplet (forced by both asymmetric pressure and shear stresses in the lubrication layer) which would significantly add to the modelling complexity, and which will not be relevant for the vertical (axisymmetric) impacts we shall consider next. We also expect these effects to be small for quasi-normal impacts (e.g. Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2019).

2.4. The lubrication-mediated droplet-bath system

We can now propose a closed system of equations for the evolution of the thin film, liquid bath and droplet. Within the bath and the droplet one may reduce the system further using the framework of Dirichlet to Neumann (DtN) maps. Consider the ‘trace’ of the potentials on the boundary of the domains in which Laplace's equation are solved:

(2.24a,b)\begin{equation} \varPhi_b(\boldsymbol{x_H},t) = \phi_b|_{z=0}, \quad \varPhi_d(\theta,\varphi,t) = \phi_d|_{r=R_0}. \end{equation}

The DtN map, denoted $D$, is the linear map that provides the normal derivative of the potential given its value on the boundary. In particular we denote

(2.25a,b)\begin{equation} (\partial_z \phi_b)|_{z=0} = D_b \varPhi_b, \quad (\partial_r \phi_d)|_{r=R_0} = D_d\varPhi_d. \end{equation}

We shall see that these maps are easily expressed in terms of eigenfunction expansions. Making use of DtN maps, we now have the evolution system

(2.26)$$\begin{gather} \partial_t \varPhi_b ={-}\frac{1}{\rho_b} p_a -g\eta_b+ \frac{\sigma_b}{\rho_b}\Delta_H \eta_b + 2\nu_b \Delta_H\varPhi_b , \end{gather}$$

(2.27)$$\begin{gather}\partial_t \eta_b = D_b \varPhi_b + 2\nu_b \Delta_H \eta_b \doteq F_b, \end{gather}$$

(2.28)$$\begin{gather}\partial_t \varPhi_d ={-}\frac{1}{\rho_d} p_a - \frac{\sigma_d}{\rho_d} \kappa_d - 2\nu_d \partial_r^2 \varPhi_d, \end{gather}$$

(2.29)$$\begin{gather}\partial_t \zeta_d = D_d \varPhi_d + 2\nu_d\left(\frac{1}{r^2}\Delta_s \zeta_d - \frac{1}{r^3} \Delta_s D_d^{{-}1} \zeta_d \right) \doteq F_d , \end{gather}$$

(2.30)$$\begin{gather}m\ddot{\boldsymbol{X}} = F ={-}\int_{\varOmega_L} p_a \, \boldsymbol{n_d} \,\text{d}a - mg \, \hat{\boldsymbol{z}}, \end{gather}$$

(2.31)$$\begin{gather}\boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{Q} = W+F_d-F_b, \quad \boldsymbol{x_H} \in \varOmega_{L}. \end{gather}$$

In the surface potential equations, the pressure $p_a$ is non-zero only in $\varOmega _L$ and $W={{\rm d}\kern0.7pt\boldsymbol {X}}/{{\rm d}t}$ is the vertical velocity of the droplet's centre of mass. In the kinematic condition for the droplet, we used the leading-order balance $\partial _t \zeta _d = D_d \varPhi _d$ to express the right-hand side locally in time. While $D_d$ was defined using the velocity potential, its inverse acts on $\zeta _d$ here. This is admissible since it is an invertible linear operator on functions with mean zero, and we may restrict the inverse map to have mean zero also. Equation (2.31), where $\boldsymbol {Q}$ is given by (2.13), is an inhomogeneous elliptic equation for the pressure with the boundary condition that $p_a=0$ on the free boundary $\partial \varOmega _L$.

3.1. Numerical implementation

We shall compute a solution approximating the system (3.3)–(3.8) by truncating the Bessel expansions at a large $N$ and integrating the resulting differential equations (3.3)–(3.5) with a fourth-order Runge–Kutta scheme. The equation for the pressure (3.8) is calculated using a trapezium rule. Since Bessel expansions are ill-conditioned on regular grids, we use an oversampled non-uniform grid in $r$ with $M$ points distributed on Chebyshev collocation points on $[0,L]$ using $\theta _j = (j-1){\rm \pi} /(M-1)$ with $j = 1\ldots M$ and $r_j = {L(1+\cos (\theta _j))}/{2}$. Hence we compute two matrices: an $M\times N$ matrix which evaluates an $N$-term Bessel series at $r_j$ and an $N\times M$ matrix corresponding to calculating the projection formula (3.7) for the Bessel coefficients of a function given at $r_j$.

3.2. Model results

We first present results of the rebound of two representative spheres of different radii and densities rebounding off a deep-water bath against a sweep of initial downward velocities. The choice of parameters is displayed in table 1 and were chosen to correspond to the data for low-velocity impacts in Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). This will permit the results of the LM model derived within this paper to be compared with the KM model (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2017) and DNS of the Navier–Stokes equations. In that study, experiments were also performed using hydrophobically coated solid spheres.

Table 1. Parameters used in simulations of a water bath and various solid spheres.

For such studies, the coefficient of restitution is usually taken as the negative ratio of the speeds between the time of impact $t_{imp}$, taken to be when the south pole of the sphere first crosses $z=0$ (correspondingly the centre of mass crosses $z=R_0$), and the time of lift-off $t_{lift}$, which is when the sphere leaves the same height on an upwards trajectory. The time spent below $z=0$ is the contact time $t_c = t_{lift}-t_{imp}$. In Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) for low impact velocities, the square coefficient of restitution $\alpha ^2$, which is the ratio of mechanical energies before and after impact (using the $z=0$ crossing as reference height), is used, as it can capture very small rebounds where $\alpha ^2<0$ where the sphere detaches but it's centre of mass does not reach the $z=R_0$ after rebound. Thus

(3.9)\begin{equation} \alpha^2 = \frac{E_{out}}{E_{in}} = \frac{W_{imp}^2}{W_{detach}^2-2g(Z-R_0)}, \end{equation}

where $W_{imp}$ and $W_{detach}$ are the vertical velocities at the moment of impact $t_{imp}$ and separation $t_{detach}$, when $\text {min}(h)>\varepsilon$ is first achieved after impact. In figure 1, $\alpha ^2$, the maximum penetration depth $\delta = -(Z_{min}-R_0)$, and the pressing time $t_p=t_{detach}-t_{imp}$ are displayed for the LM and KM models and DNS. In Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) $t_p$ is favoured over $t_c$ to capture the rebounds which do not return to $z=0$.

Figure 1. Comparison of $\alpha ^2$, $\delta$ and $t_p$ for two parameter regimes, with LM ($\circ$), compared with DNS ($\Diamond$) and KM ($*$) from Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021).

Figure 1 demonstrates that the agreement is good between the LM model against both the results from the KM and DNS, with general trends being upheld. We note here that the latter two methods also involve certain approximations. The KM model assumes no lubrication layer (the method matches the position and velocities of the free surface to those of the sphere) and uses a tangential contact between the sphere and the bath throughout the impact. In addition, for the KM model we are comparing with, the evolution equations for the bath are linear and with the same quasi-potential approximation that we use here. The DNS results in Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) were performed on the open source Gerris where the solid drop is modelled by using artificially high viscosity and surface tension coefficients. Gerris uses a volume of fluid formulation and interfaces are approximated using continuous density changes (see Popinet Reference Popinet2009).

We conclude by providing numerical results for a configuration which can be compared with the KM model, DNS and also experimental values, as shown in table 2. The experimental configuration used a rectangular container while our simulation is cylindrical. In table 2 we varied the radius of the bath over comparable values to their length and width dimensions, resulting in a range of values for our simulations.

Table 2. Results for a falling sphere of radius $R_0= 0.83\times 10^{-4}$ m, density $\rho = 1200$ kg m$^{-3}$ and initial velocity $W_0 = -0.35\,{\rm m}\,{\rm s}^{-1}$ onto a deep-water bath. The range in the LM data arises from varying bath radius $L$.

The results for the LM model restitution coefficient $\alpha$ is at the lower end of the experimental range, the depth of impact is somewhat below the experimental range and the impact time is well within the experimental range. In addition to modelling assumptions, one potential source of error is that the physics of the hydrophobic coating used in experiments is not considered in any of the models.

Turning to impact details for this case, figures 2(a) and 3(b), show three distinct phases of rebound. First, a rapid expansion of the air layer, driven by the inertia of the impacting drop; second, a slow contraction of the layer; third, a rapid suction event due to the retraction of the free surface. The expansion of the region lasts $O(1\,\text {ms})$, the contraction lasts $O(10\,\text {ms})$, and the suction event lasts $O(1\,\text {ms})$ and precedes separation. A highlight of the LM model is the ability to capture and quantify the behaviour within the air layer during the rebound. In particular, we can study the evolution of the pressure profile and corresponding height layer within the system, quantities which are challenging to extract even in the DNS setting. In what follows, we have non-dimensionalised all lengths with the radius of the sphere (e.g. $r^*=r/R_0$, and area $A^* = A/R_0^2$), the pressure with the surface tension pressure jump across a spherical interface $P^* = P/(2\sigma R_0^{-1})$, and the force with $F^* = F/(2\sigma R_0)$. Throughout the impact, the height of the air layer retains thickness $O(1\,\mathrm {\mu }{\rm m})$, as expected from the literature (Tang et al. Reference Tang, Saha, Law and Sun2019). However, there is some variation as shown in figure 2(a), with a thickness of 0.5–2 $\mathrm {\mu }$m in the bulk of the layer and a narrowing on its boundary to 0.15–0.5 $\mathrm {\mu }$m. In effect, for the bulk of the impact, the air layer forms a quasi-static cushion constricted on its circular boundary. This narrowing effect is enhanced when the lubrication region boundary ‘turns around’ and immediately prior to detachment, where the lubrication layer has shrunk to a small disc at the base of the sphere. Turning to the profile of the pressure within the air layer, shown in figure 3(b), we find a range of 0.9–1.5 $2\sigma R_0^{-1}$ over the vast majority of the region. The maximum pressure is felt at impact and is 7.6 $2\sigma R_0^{-1}$ (indicating an inertially driven pressure) and the minimal pressure being $-$3.8 $2\sigma R_0^{-1}$ at liftoff. A pressure spike occurs where the layer thins at its edge, aligning with results found in Hendrix et al. (Reference Hendrix, Bouwhuis, van der Meer, Lohse and Snoeijer2016). Figure 3(a) captures the pressure profiles during expansion and contraction. Small oscillations observed in the pressure profile at two intermediate times, we believe, are numerical artefacts, as they are lessened by increased mesh resolution. Lesser refined meshes had larger oscillations, however, with minimal effect on the overall dynamics (restitution coefficient, impact times, generated waves, etc.) of the system.

Figure 2. (a) Thickness of the air layer, $h(r^*,t)$. (b) Free surface of the bath, $\eta ^*_b(r^*,t)$, shows the sinking depth of the sphere relative to the resultant waves on the free surface.

Figure 3. (a) Lubrication-layer pressure profiles $P^*(r^*,t)$ at selected times. Dark and light curves correspond to expanding and contracting layers, respectively. (b) Colourmap showing near uniform pressure in the interior and its rapid expansion and slow contraction.

Figure 4 shows a temporal slice of pressure, under the south pole of the sphere $P(0,t)$, and the total pressure force acting on the sphere. The pressure clearly has fine-scale behaviour at the initial impact and detachment but these features do not substantially affect the overall force on the system. Most of the work on the drop is due to the quasi-uniform pressure at intermediate times. Figure 2(b) shows the full wavefield during impact.

Figure 4. Comparison between the pressure force $F^*$ felt by the droplet, the pressure $P^*$ along $r=0$ during impact, and the pressed area $A_p^*$. Note initial (inertial) and final (suction) spikes. The consequences of the force on the wider system is well tracked by the lubrication-layer area.

All computations presented were performed in Matlab with ranges $N=2^{11}-2^{12}$, $M=2^8-2^9$. For lower impact velocities this was an over-resolution, whereas the finer grids were needed for stronger impacts. The domain length $L$ varied from $12R_0$ to $16R_0$, the latter ensuring the impact is unaffected by waves reflecting from the boundary, even for longer impacts. The longest simulations took $O(10\,\text {h})$ on a standard personal computer running Matlab, although the use of an adaptive time step would probably reduce this considerably. The only free parameter in the system is $\varepsilon$, the cutoff for lubrication, and this was set to ${\rm 40}\,\mathrm {\mu }{\rm m}$, which corresponds to approximately 5 % the radius of the sphere $R_0$, though values between 10 % and 1 % were also tested and the model displayed insensitivity within this range of $\varepsilon$. The fast decay in pressure in (3.8) as $h$ increases renders the solution insensitive to larger values of the cutoff $\varepsilon$.