1. Introduction
The dynamics of a falling liquid drop as it impacts on a substrate depend on the initial conditions and physical properties of the liquid, the surrounding media and the substrate. Depending on the relevant physical parameters, the liquid drop may splash, spread or bounce upon impact with a solid substrate. The diversity of patterns produced by liquid drops impacting on a solid substrate in the presence of ambient air were first captured by Worthington (Reference Worthington1877), even before the invention of flash photography, through observations of patterns under the light of a spark produced by a break in an electric circuit when the drop hits the substrate. Worthington's drawings of these patterns sparked many investigations into the physical phenomena associated with drop impact, many of which are summarized in review articles (Rein Reference Rein1993; Yarin Reference Yarin2006; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016).
Many physical parameters are known to affect the dynamics of a drop upon impact with a substrate (Rioboo, Marengo & Tropea Reference Rioboo, Marengo and Tropea2002). These include the initial size and velocity of the drop (Pasandideh-Fard et al. Reference Pasandideh-Fard, Qiao, Chandra and Mostaghimi1996; Riboux & Gordillo Reference Riboux and Gordillo2014), viscosity of the liquid (Pasandideh-Fard et al. Reference Pasandideh-Fard, Qiao, Chandra and Mostaghimi1996; Stevens, Latka & Nagel Reference Stevens, Latka and Nagel2014), pressure and viscosity of the surrounding medium (Xu, Zhang & Nagel Reference Xu, Zhang and Nagel2005; Sprittles Reference Sprittles2015), interfacial tension between the liquid and its surrounding medium (Pasandideh-Fard et al. Reference Pasandideh-Fard, Qiao, Chandra and Mostaghimi1996; Rioboo et al. Reference Rioboo, Bauthier, Conti, Voue and De Coninck2003), and the physical properties of the substrate (Rein Reference Rein1993; Range & Feuillebois Reference Range and Feuillebois1998; Yarin Reference Yarin2006; Latka et al. Reference Latka, Strandburg-Peshkin, Driscoll, Stevens and Nagel2012; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016). Experiments by Xu et al. (Reference Xu, Zhang and Nagel2005) show that in the impact of a liquid drop on a solid substrate, ambient air pressure determines whether the drop eventually splashes, characterized by the ejection of a number of small drops into the air, or forms a thin film and spreads smoothly onto the surface (Driscoll, Stevens & Nagel Reference Driscoll, Stevens and Nagel2010). This suggests a role of the surrounding gaseous medium in determining the impact dynamics even before the drop makes contact with the solid substrate, motivating investigation into the role of an intervening gaseous layer in determining dynamics of a liquid drop as it approaches a solid substrate, prior to contact of the drop with the substrate.
Theoretical work on the role of an intervening layer of gas in the approach of a liquid drop towards a solid surface has been carried out in simplified two-dimensional (Mandre, Mani & Brenner Reference Mandre, Mani and Brenner2009; Mani, Mandre & Brenner Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012; Moore, Ockendon & Oliver Reference Moore, Ockendon and Oliver2013) and axisymmetric (Hicks & Purvis Reference Hicks and Purvis2010) geometries. Some of these theoretical investigations (Mandre et al. Reference Mandre, Mani and Brenner2009; Hicks & Purvis Reference Hicks and Purvis2010; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012) consider an inviscid liquid drop falling towards a solid surface in the presence of an ideal gas, and predict that the drop deforms due to a buildup of pressure in the gas trapped underneath the drop, prior to its contact with the solid substrate. Since these studies assume the fluid flow in the drop is inviscid, the velocity field is given by a potential flow. This substantially simplifies the analysis, allowing the dynamics to be modelled using values of relevant field variables exclusively at the liquid–gas interface using a lubrication approximation, without explicitly computing the fluid flow in the bulk of the drop. A schematic of the deformation of the drop is shown in figures 1(a) and 1(b).
Schematic of the physical and computational model. (a) Relative locations of the falling drop, intervening gas layer and solid substrate. (b) A schematic zoomed into the region of interest, showing an exaggerated view of a deformed liquid–gas interface (green), and a schematic of the region where viscous effects are important (yellow). (c) Control volume and grid for discretization constructed in the fluid domain.
The theoretical predictions of drops deforming on a layer of air between the liquid drop and solid substrate, prior to contact between the drop and the substrate, are consistent with previous experimental observations (Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005). Measurements of the air layer thickness underneath the impacting drop have been performed (Driscoll & Nagel Reference Driscoll and Nagel2011; Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; van der Veen et al. Reference van der Veen, Tran, Lohse and Sun2012). Subsequent experimental work by Kolinski, Mahadevan & Rubinstein (Reference Kolinski, Mahadevan and Rubinstein2014b) shows the effect of liquid viscosity on the evolution of the liquid–air interface even before contact with the substrate, which is not captured by the potential flow description of the evolution of liquid dynamics in previous theoretical work (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012).
In this work we are concerned with the dynamics of a liquid drop as it approaches a rigid, non-porous, solid surface, in the presence of an ambient gas. We modify the physical model of Mandre, Mani and Brenner by incorporating the viscosity of the liquid. Due to the additional viscous term in the governing equations for the liquid, the coupled interaction of the liquid and gas becomes theoretically intractable, and it is necessary to solve numerically for the fluid flow in the bulk of the drop. There are a variety of previous computational studies of drop impact, such as looking at the dynamics of spreading (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010), splashing (Josserand & Zaleski Reference Josserand and Zaleski2003; Boelens & de Pablo Reference Boelens and de Pablo2018) or rebound from superhydrophobic surfaces (Renardy et al. Reference Renardy, Popinet, Duchemin, Renardy, Zaleski, Josserand, Drumright-Clarke, Richard, Clanet and Quéré2003). Duchemin & Josserand (Reference Duchemin and Josserand2011) generalized the potential flow model of Mandre et al. (Reference Mandre, Mani and Brenner2009) to work in axisymmetric coordinates, where the liquid–gas interface is represented as an arbitrary curve, allowing a full spherical drop to be represented. With this model, they were able to predict the emergence of a thin jet skating above the gas layer, and they found that without surface tension a finite-time singularity forms where the liquid–gas interface touches the substrate. More recently, Duchemin & Josserand (Reference Duchemin and Josserand2020) solved the lubrication equations in the thin gas layer; this work describes the drainage of the gas during impact of a drop onto a solid substrate or a liquid film.
Philippi, Lagrée & Antkowiak (Reference Philippi, Lagrée and Antkowiak2016) and Jian et al. (Reference Jian, Josserand, Popinet, Ray and Zaleski2018) examined the early stages of drop impact, making use of the Gerris/Basilisk flow solvers (Popinet Reference Popinet2003, Reference Popinet2009; Lagrée, Staron & Popinet Reference Lagrée, Staron and Popinet2011; Popinet Reference Popinet2015) for tracking the liquid–gas interface. This simulation framework provides greater flexibility allowing the dynamics to be tracked over a longer time. In comparison, our reduced description uses a simpler description of the initial dynamics, and is well suited to handle the disparate length scales between the gas layer and the drop. We solve the incompressible Navier–Stokes equations in the drop, using a modern implementation of Chorin's projection method (Chorin Reference Chorin1967, Reference Chorin1968) that incorporates improvements from Almgren, Bell, Collela and coworkers (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989; Colella Reference Colella1990; Almgren, Bell & Szymczak Reference Almgren, Bell and Szymczak1996). Our numerical model allows us to resolve the flow field and pressure in the drop, and examine the effects of many different physical parameters in ways (e.g. viscosity, surface tension) that would be difficult to do in an experiment. We validate our model using theoretical results (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012). Using our model, we are able to recapitulate the square-root scaling with liquid viscosity of the lift-off time observed by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b). However, our results show that the precise relationship between lift-off time and liquid viscosity is more complicated. We explore the dependence of lift-off time on surface tension, initial drop velocity and drop radius.
In the following sections we describe the physical model and parameter regime, followed by approximations made in the simulation domain. We then describe the results from our simulations, along with comparisons with previous theoretical and experimental studies. Our numerical results are consistent with theoretical calculations (Mandre et al. Reference Mandre, Mani and Brenner2009) as well as experimental results (Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014b). Our results predict the effect of viscosity and interfacial tension on the dynamics of drop impact.
2. Model
In this section we first explain the physical set-up and specify the parameter regime considered in the rest of this work. Based on the physical set-up and parameter regime, we then describe the mathematical model. The mathematical model largely derives from the previous theoretical work of Mandre et al. (Reference Mandre, Mani and Brenner2009), Mani et al. (Reference Mani, Mandre and Brenner2010) and Mandre & Brenner (Reference Mandre and Brenner2012). We specify the predictions made by this model, experimentally observed deviations from these predictions, and then the mathematical model used in our work. We then use this mathematical model to derive appropriate boundary conditions for the flow in the drop.
2.1. Physical problem and parameter regime
The physical set-up is a drop of liquid falling towards a flat, solid surface in the presence of a surrounding gas. As the liquid drop falls towards the solid surface, it interacts with the surrounding gas. We are interested in modelling and simulating the coupled dynamics of the liquid and the gas before the liquid makes contact with the solid surface.
Figure 1(a) shows a schematic of the physical problem, indicating the relative locations of the liquid drop, the surrounding gas and the solid surface. The drop is initially spherical, with radius 
$R$, falling towards a horizontal solid surface at uniform vertical velocity 
$V_0$. The gas is initially at uniform pressure 
$P_0$. We specify the values of initial conditions and physical properties considered for the liquid and air in table 1. We use the subscripts 
$l$ and 
$g$ to denote the liquid drop and the gas, respectively, and 
$0$ to denote the initial conditions.
Relevant initial conditions and physical parameters, and their approximate values, which informs the mathematical model. The subscripts 
$l$ and 
$g$ denote the liquid drop and the gas, respectively, and 
$0$ denotes initial conditions.
In the parameter regime specified in table 1 potential flow theory argues that the relevant Reynolds number, based on the initial velocity 
$V_0$ and drop radius 
$R$, is 
${O}(10^{2})$, allowing for an inviscid consideration of the liquid. We refer to the mathematical model as the potential flow model, and the theoretical predictions arising from this mathematical model as the potential flow theory.
The potential flow theory predicts that, as the drop falls towards the solid surface, a buildup of air pressure underneath the drop causes the drop to deform in the middle, developing a dimple. The drop then spreads over a thin layer of gas that separates it from the substrate. The potential flow theory predicts scaling laws for the height of the gas layer. Subsequent experimental observations of the gas underneath the drop (Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; van der Veen et al. Reference van der Veen, Tran, Lohse and Sun2012) show the existence of this dimple and the shape of the drop profile, similar to the predictions of the potential flow theory.
Experimental observations made by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) in part of this parameter regime show that the falling drop first develops a dimple in the middle of the drop, as schematized by the liquid–gas interface shown in green in figure 1(b). This is consistent with the potential flow theory (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012). Subsequently in time, and prior to contact of the liquid–gas interface with the solid substrate, the shape of the interface depends on the kinematic viscosity of the liquid. Specifically, there is a critical time 
$\tau _c$ at which the minimum height of the drop from the substrate stops decreasing, and the leading edge of the drop begins to move away from the surface. The time 
$\tau _c$ is measured relative to the time when the drop would have made contact with the substrate if it were not deforming due to interactions with an intervening gas. Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) observe that 
$\tau _c \propto {\nu _l}^{1/2}$.
As a consequence of the discrepancy between the observations of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b), who definitively show the effect of liquid viscosity on the shape of the liquid–gas interface well before contact between the liquid and the substrate, and the potential flow theory, we are specifically interested in capturing the role of liquid viscosity, on these coupled dynamics. Based on the observations of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b), we modify the potential flow model to consider a liquid described by the full Navier–Stokes equations, rather than an inviscid approximation.
2.2. Mathematical model
The dynamics of an initially spherical drop are naturally described in terms of three-dimensional spherical polar coordinates, and those of the gas layer in terms of three-dimensional cylindrical polar coordinates. In this coordinate system, based on experimental observations of air layer profiles under an impacting drop (Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; van der Veen et al. Reference van der Veen, Tran, Lohse and Sun2012; Kolinski, Mahadevan & Rubinstein Reference Kolinski, Mahadevan and Rubinstein2014a; Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014b, Reference Kolinski, Kaviani, Hade and Rubinstein2019), we first use the simplifying approximation of axisymmetric flow in the drop and the gas, with the axis of symmetry being perpendicular to the substrate, through the centre of the undeformed drop.
As in the potential flow model (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012), we further simplify our consideration by approximating the drop as initially being an infinite cylinder, with the long axis perpendicular to the substrate, rather than spherical. With these assumptions and approximations, we can describe the dynamics of the liquid and gas using two-dimensional Cartesian coordinates 
$(x,y)$. The initial height 
$h(x,t=T_0)$ of the liquid–gas interface is described by a parabolic profile as
\begin{equation} h(x,t=T_0) = H_0 + \frac{x^2}{2R}, \end{equation}
where 
$H_0$ is the initial height of the centre of the interface, i.e. at 
$x=0$.
We model the flow in the liquid using the incompressible Navier–Stokes equations. The continuity equation for the liquid is
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_l = 0, \end{equation}and the momentum equation for the liquid is
\begin{equation} \boldsymbol{u}_{l,t} + \boldsymbol{u}_l \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}_l ={-}\frac{\boldsymbol{\nabla}p_l}{\rho_l} + \nu_l \nabla^2 \boldsymbol{u}_l, \end{equation}
where 
$\boldsymbol {u}_l = (u_x, u_y)$, 
$p_l$ and 
$\rho _l$ are the velocity, pressure and density of the liquid, respectively.
Based on the potential flow theory (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012) and subsequent experimental observations (de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; van der Veen et al. Reference van der Veen, Tran, Lohse and Sun2012; Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014b), the gas layer is slender, with the horizontal length scale 
$L_x$ being much greater than the vertical length scale 
$L_y$, allowing for an approximation of one-dimensional flow in the gas. The height 
$h=h(x,t)$ of the gas layer is thus described by the compressible lubrication equation (Taylor & Saffman Reference Taylor and Saffman1957),
\begin{equation} 12 \mu_g (\rho_g h)_t = (\rho_g h^3 p_{g,x})_x, \end{equation}and the equation of state is described by an adiabatic assumption,
\begin{equation} \frac{p_g}{P_0} = \left(\frac{\rho_g}{\rho_0}\right)^{\gamma}. \end{equation}The coupling of flows in the liquid and gas at their interface is done using the Laplace condition for pressure,
\begin{equation} p_l = p_g + \sigma h_{xx}, \end{equation}and equality of shear stress in the liquid and gas at the liquid–gas interface is
\begin{equation} \tau_l(x,h) = \tau_g(x,h). \end{equation}
Far enough away from the deforming effects of the intervening gas layer on the liquid drop, the pressure in the gas is the ambient pressure, 
$\lim _{x\rightarrow \pm \infty } p_g(x,t) = P_0$.
Given that the initial and boundary conditions are symmetric about 
$x=0$, we modify our mathematical model to incorporate symmetry boundary conditions about 
$x=0$. We do so by applying the boundary conditions 
$(u,v_x)=(0,0)$ in the liquid domain and 
$p_{g,x}=0$ and 
$h_x=0$ in the gas layer.
2.3. Liquid domain and boundary conditions
Having specified the physical problem in § 2.1 and the mathematical model, with initial and boundary conditions, in § 2.2, we now turn to a description of the domain and boundary conditions that specify the flow in the liquid in a computationally tractable manner. The assumptions that lead to numerical approximations in this section are specified in §§ 2.1 and 2.2. We detail them here as appropriate.
2.3.1. Rectangular domain with mass flux
Following from the previous section, the liquid drop is described by a two-dimensional flow field. As the drop interacts with the gas layer, its shape changes. The drop begins to deform close to the middle of the interface, at 
$x=0$. With this observation, we choose a region of interest, as shown in figure 1(b), that stretches outwards from 
$x=0$ and upwards from 
$y=h(x,t)$.
We make an additional approximation that the flow in the liquid at 
$y=h(x,t)$ can be approximated by the flow at 
$y=0$. We can remove the spatial dependence because the liquid–gas interface is slender. We simulate the control volume with and without a temporal dependence 
$y=h(t)$, and do not observe a significant change in results. With these approximations, we are able to choose a rectangular, inertial control volume, stretching outwards from 
$x=0$ and upwards from 
$y=0$, whose size is chosen based on the relevant physical scales (§ 2.4). With the use of symmetry in the domain, as described in § 2.2, we need to model only half of the control volume. We choose to model the right half, leading to the computational domain shown in figure 1(c). In subsequent sections, we describe the boundary conditions at each of the boundaries in the computational domain.
2.3.2. Treatment of the viscous term and bottom boundary condition
Initially, the drop is moving at uniform vertical velocity, meaning that all spatial velocity gradients are zero, and 
$\nabla ^2 \boldsymbol {u}_l=0$. Deformation of the drop begins at 
$x=0$. Our model consequently assumes that the shear stresses are small away from this region where the drop has deformed. In a two-dimensional flow vorticity is introduced only through shear stresses at the boundary. The vorticity therefore spreads into the interior of the drop from the deformed interface. We thus model the drop as having a region close to 
$(x,y)=(0,0)$, where viscosity is important, while far away from this region, 
$\nu _l\nabla ^2\boldsymbol {u}_l$ makes a negligible contribution to (2.3). A schematic of the region where viscosity is important in the deformed drop is shown by yellow in figures 1(b) and 1(c). Consequently, we model the domain as having one boundary where viscous effects are important, which is the bottom boundary.
At the bottom boundary, we apply the boundary conditions specified in (2.6) and (2.7). Equation (2.7) is enforced through a condition on the gradient of the horizontal velocity,
\begin{equation} u_{l,y}(x,0)=\frac{\mu_g}{\mu_l}u_{g,y}(x,0). \end{equation}2.3.3. Boundaries in the interior of the drop: top and right boundary conditions
According to our approximation of a rectangular domain, the bottom boundary of the liquid interacts with the gas layer. The other boundaries are chosen to be in the interior of the drop, away from the effect of vorticity from the bottom boundary, and can be treated as inviscid. With this approximation, at these boundaries, denoted as interior boundaries, we simplify the momentum equation for the liquid, (2.3), by noting that the viscous term, 
$\nu _l\nabla ^2\boldsymbol {u}_l$, is relatively small.
Additionally, the velocity gradients in the liquid are initially zero. Away from the effects of the gas layer, the nonlinear advective term in the momentum equation, 
${\boldsymbol u}_l \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u}_l$ in (2.3), will also be small. With these approximations, we obtain a simplified momentum equation for the liquid at these boundaries,
\begin{equation} \boldsymbol{u}_{l,t} ={-}\frac{\boldsymbol{\nabla}p_l}{\rho_l}. \end{equation}
We use this simplified momentum equation to obtain velocity boundary conditions at the interior boundaries. Taking the divergence of (2.9) and noting that the flow in the liquid is incompressible, we obtain 
$\nabla ^2 p_l = 0$. Knowing the pressure at the bottom boundary, and approximating the liquid as a semi-infinite domain, we get an analytical expression for the pressure at the interior boundaries,
\begin{equation} p_l(x,y) = \frac{1}{\rm \pi}\int_{-\infty}^{\infty}p_l(s,0) \frac{y}{(x-s)^2+y^2} \, \mathrm{d}s. \end{equation}
Assuming that the domain is larger than the region where the pressure from the gas is non-zero, we numerically integrate over a finite length of the domain boundary to get the pressure at a particular 
$(x,y)$ location along the interior boundaries. Given the pressure at the boundary, we integrate (2.9) in time to obtain velocity boundary conditions at the interior boundaries. The size of the region where vorticity is significant determines the size of the domain, and the locations of the interior boundaries, which are calibrated in simulation.
2.4. Choice of domain
We simulate the coupled dynamics of the liquid and gas over the time interval 
$0\le t \le t^{end}$, for a simulation domain of 
$0 \le x \le L$ and 
$0 \le y \le \beta L$, where 
$\beta$ is the aspect ratio, for the liquid, and a domain of 
$0 \le x \le L$ for the gas. To set the size of the domain and duration of the simulation, we consider the initial dynamics as the drop falls toward the surface, and compare to the theoretical results of the potential flow theory (Mandre et al. Reference Mandre, Mani and Brenner2009). We use the scaling results of the potential flow theory for the centre of the gas layer, 
$h(x=0,t)$. The scaling analyses predict the height 
$H^*$ at which the pressure buildup in the gas layer causes the drop to undergo significant deformation and develop a stagnation point at 
$x=0$, forming a dimple.
According to potential flow theory, the compressibility of the flow in the gas is determined by the parameter
\begin{equation} \epsilon \equiv \frac{P_0 R \, St^{4/3}}{\mu_g V}, \end{equation}
which is the ratio of the initial gas pressure to the pressure that would have built up if the gas were treated as incompressible, where 
$St \equiv {\mu _g}/{\rho _l V R}$ (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012; Langley, Li & Thoroddsen Reference Langley, Li and Thoroddsen2017).
The results of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b), showing the relationship between 
$\nu _l$ and the drop profile, correspond to 
$\epsilon ^{-1} < 1$, meaning that the flow in the gas can be considered to be incompressible. While our simulations model the full incompressible flow in the gas, we use this approximation to set the initial conditions. In this regime we have the estimate
\begin{equation} H^* = R \, St^{2/3}. \end{equation}
This scaling result is strictly valid for inviscid flow, 
$\nu _l = 0$, 
$Re \rightarrow {}\infty$, and no surface tension at the liquid–air interface, 
$\sigma = 0$, 
$We \rightarrow {}\infty$. While our simulations incorporate liquid viscosity and surface tension, (2.12) still provides a good estimate for the height at which the stagnation point forms, and is a useful measure of the physical scales of interest. We therefore set the initial height to be a multiple of 
$H^*$, so that
\begin{equation} H_0 = \tilde{H}_0 H^* = \tilde{H}_0 R \, St^{2/3} \end{equation}
for a dimensionless constant 
$\tilde {H}_0$. Similarly, we set 
$\tilde {t}_{end}= \tilde {t}_{end} R \, St^{2/3}/V$, and based on the parabolic shape of the initial profile we choose 
$L= \tilde {L} R \, St^{1/3}$, where 
$\tilde {L}$ and 
$\tilde {t}_{end}$ are dimensionless constants.
3. Numerical implementation
In this section we describe the numerical implementation of the mathematical model described in § 2.2. Throughout the liquid domain, we compute and record the liquid velocity 
$\boldsymbol {u}_l$ and pressure 
$p_l$. In the gas layer we compute and track the height 
$h$ and pressure 
$p_g$. A reader may skip the remainder of this section without loss of continuity.
The flow-field variables 
$\phi$ are computed at discretized timesteps, with simulation time 
$t=n {\rm \Delta} t$, where 
${\rm \Delta} t$ is the timestep discretization used in the simulation and 
$n$ is an integer counter. Given values of field variables 
$\phi ^{(n)}=\phi (t=n {\rm \Delta} t)$, the goal of the simulation is to compute the field variables 
$\phi ^{(n+1)}$. This section describes the numerical method for doing so, for the liquid flow-field variables, 
$\boldsymbol {u}_l$ and 
$p_l$, and gas flow-field variables, 
$h$ and 
$p_g$.
3.1. Projection method
The core of the algorithm is based on Chorin's projection method (Chorin Reference Chorin1967, Reference Chorin1968), with further developments in the computational techniques behind the fluid simulation method (Bell et al. Reference Bell, Colella and Glaz1989; Almgren et al. Reference Almgren, Bell, Colella, Howell and Welcome1998). In particular, the work of Almgren et al. (Reference Almgren, Bell and Szymczak1996) introduces the approximate projection method discretized using the finite-element method (FEM). The numerical methods are described in detail by Sussman et al. (Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999), Yu, Sakai & Sethian (Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007) and Rycroft et al. (Reference Rycroft, Wu, Yu and Kamrin2020). Here, we sketch the main features of these methods.
We solve for 
$\boldsymbol {u}_l$ and pressure 
$p_l$ in the liquid domain by solving the momentum equation for the liquid, (2.3), subject to incompressibility of the liquid, as specified in (2.2). In describing the method for the solution of the field variables in the liquid domain, for the remainder of this section, we drop the subscript 
$l$.
The field variables are defined on a two-dimensional grid, discretized into rectangular cells of size 
${\rm \Delta} x$ by 
${\rm \Delta} y$, as shown in figure 2. The velocities 
$\boldsymbol {u}^n$ are located in the cell centres, and the pressures 
$p^n$ are located at the cell corners. To advance from step 
$n$ to 
$n+1$, we first compute an intermediate velocity 
$\boldsymbol {u}^*$ according to
\begin{equation} \frac{{\boldsymbol u}^{*}-{\boldsymbol u}^n}{{\rm \Delta} t} ={-} [{\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u}]^{n+1/2} + \frac{\nu}{2} (\nabla^2 {\boldsymbol u}^n + \nabla^2 {\boldsymbol u}^* ). \end{equation}
Here, the viscous stress terms 
$\nu \nabla ^2 \boldsymbol {u}$ and 
$\nu \nabla ^2 \boldsymbol {u}^*$ are evaluated using a standard five-point finite-difference stencil. Due to the presence of the 
${\boldsymbol u}^*$ on the right-hand side of (3.1), it must be solved implicitly. We use a multigrid method to solve the resulting linear system.
Discretization of the field variables. (a) The two-dimensional computational grid for the liquid. (b) The one-dimensional computational grid for the gas.
The advective term 
$[\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}]^{n+1/2}$ is evaluated at the half-timestep 
$n+1/2$ using the second-order Godunov upwinding scheme of Colella (Reference Colella1990). This is performed by extrapolating the velocity in each cell to midpoints of the four edges using a first-order Taylor expansion. Doing this requires that the gradient of the velocity is first evaluated at each cell centre, which is done using the fourth-order monotonicity-limited scheme of Colella (Reference Colella1985). This scheme uses a five-point stencil in each coordinate, therefore requiring information from two grid points in each orthogonal direction.
After this procedure, each edge has two velocities from the two adjacent cells. A Godunov upwinding procedure is then used to select one based on the direction of the velocity. An intermediate marker-and-cell (MAC) projection is then used to adjust the velocities to satisfy a discrete zero divergence criterion, so that the net mass flow out of each cell is zero (Bell, Howell & Colella Reference Bell, Colella and Howell1991; Sussman et al. Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999). After this, the advective term can be accurately evaluated using centred finite differences of the upwinded edge-based velocities (Sussman et al. Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999; Yu et al. Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007; Rycroft et al. Reference Rycroft, Wu, Yu and Kamrin2020).
The velocity at step 
$n+1$ can be calculated from the intermediate velocity by evaluating the pressure gradient term,
\begin{equation} \frac{{\boldsymbol u}^{n+1}-{\boldsymbol u}^*}{{\rm \Delta} t} ={-}\frac{1}{\rho}{\boldsymbol{\nabla} p^{n+1}}, \end{equation}
but this requires knowledge of the pressure field 
$p^{n+1}$, and there is no explicit update equation for this field in the incompressible limit. To proceed, we take divergence of (3.2) and enforce that 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}^{n+1}=0$, according to (2.2). Hence,
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[{\rm \Delta} t \left\{ \frac{1}{\rho} \boldsymbol{\nabla} p^{n+1} \right\}\right] = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol u}^*, \end{equation}
which is a Poisson equation for the pressure that can be solved. Once 
$p^{n+1}$ is known, (3.2) can be used to evaluate 
$\boldsymbol {u}^{n+1}$ and complete the timestep from 
$n$ to 
$n+1$. We solve (3.3) using the FEM discretization of Almgren et al. (Reference Almgren, Bell and Szymczak1996). Here, each pressure value at a cell corner has a corresponding bilinear hat function over the four neighbouring grid cells (Yu et al. Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007; Rycroft et al. Reference Rycroft, Wu, Yu and Kamrin2020). Once the pressure is computed, the gradient in (3.2) is evaluated using centred finite differences of the pressure field values.
3.2. Implementation of boundary conditions: ghost layers
The boundary conditions are implemented by means of ghost layers. Around the boundaries of the two-dimensional liquid domain, there are two layers of additional grid points. Two layers are required in order to evaluate the fourth-order monotonicity-limited derivatives arising from the advective term. The boundary conditions are specified in these layers of grid points. These conditions are then used to compute the necessary gradients in (3.1).
The velocity boundary conditions at the top and right boundaries of the liquid domain, corresponding to the interior of the drop, and described in § 2.3.3, are obtained through substituting the expression for pressure given by (2.10) in the expression for velocity at these boundaries, (2.9), and then using Simpson's rule to numerically integrate the pressure along these boundaries.
3.3. Solution of the gas layer equations
The field variables 
$h$ and 
$p_g$ for the gas layer are discretized in accordance with the discretization for the field variables for the liquid. Specifically, the pressure field 
$p_g$ is stored at grid points in the same locations along the horizontal axis as the pressure field in the liquid, as shown in figure 2. The height 
$h$ is stored at grid points in the same locations along the horizontal axis as the velocity field 
$\boldsymbol {u}$ in the liquid.
As shown in Appendix A.1, substituting the equation of state (2.5) into the lubrication equation (2.4) yields
\begin{equation} 12 \frac{\mu}{\gamma} h p_{t} + 12 \mu h_t p = \frac{1}{\gamma} h^3 p_{x} p_{x} + 3 h^2 h_x p p_{x} +h^3 p p_{xx}, \end{equation}
where the 
$g$ subscript on the pressure has been dropped. We solve for this equation as follows:
(i)
$p$ is known at the timestep $t^{(n)} = n {\rm \Delta} t$;(ii)
$h$, $h_t$ and, therefore, $h_x$ are known at timestep $t^{(n)} = n {\rm \Delta} t$;(iii) solve for
$p$ at timestep $t^{(n+1)} = (n+1) {\rm \Delta} t$;(iv) use
$p$ as a boundary condition for the flow in the liquid, to obtain $v=h_t$ at timestep $t^{(n+1)}$.
To numerically solve (3.4), we discretize the spatial derivatives using second-order centred finite differences. To avoid a timestep restriction, we make use of the backward Euler method for discretizing the temporal derivative. Hence, we obtain the discretized form
\begin{align} &12 \frac{\mu}{\gamma} \bar{h} \left(\frac{p^{n+1}_i - p^n_i}{{\rm \Delta} t} \right) + 12 \mu \bar{h}_t p^{n+1}_i= \frac{1}{\gamma} \bar{h}^3 \left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2 {\rm \Delta} x} \right)^2 \nonumber\\ &\quad + 3 \bar{h}^2 \bar{h}_x p^{n+1}_i \left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2 {\rm \Delta} x} \right) + \bar{h}^3 p^{n+1}_i \left(\frac{p^{n+1}_{i+1}-2 p^{n+1}_{i} + p^{n+1}_{i-1}}{{\rm \Delta} x^2} \right), \end{align}
that can be solved for the gas pressure 
$p^{n+1}_i$. In (3.5) the terms 
$\bar {h}$, 
$\bar {h}_t$ and 
$\bar {h}_x$ must be evaluated at the location of 
$p^{n+1}_i$. Since the height field is staggered with respect to the pressure field, this is done via linear interpolation and centred differencing, so that
\begin{equation} \bar{h} = \frac{h_i+h_{i+1}}{2}, \quad \bar{h}_t = \frac{h_{t,i} + h_{t,i+1}}{2}, \quad \bar{h}_x = \frac{h_{i+1}-h_i}{{\rm \Delta} x}. \end{equation}Due to the products of pressure terms on the right-hand side of (3.5), it is a nonlinear system of equations. We solve this using the Newton–Raphson method as described in Appendix A.2.
3.4. Code implementation and parameter choices
The simulations are performed using a custom code written in C++, using the OpenMP library for multithreading. The code makes use of the templated geometric multigrid (TGMG) library for solving the linear systems that arise when solving for the fluid. Each timestep requires four linear system solves: (1) to apply the MAC projection, (2) to apply the approximate FEM projection, and (3,4) to solve for the 
$x$ and 
$y$ velocity components when treating the viscous stress term implicitly. The code accepts a text configuration file as input, which sets all of the physical parameters, and describes the computational domain. The code and sample text configuration files are available on GitHub – see the data availability statement. The simulations are performed on a 
$M \times N$ grid for the liquid domain, and grid of length 
$M$ for the gas layer. We choose the grid spacings to be equal, so that 
${\rm \Delta} x={\rm \Delta} y$. This implies that the aspect ratio 
$\beta$ is equal to 
$N/M$.
Since the second derivative terms in both the liquid domain and the gas layer are handled implicitly, and the surface tension term in (2.6) is small, there is no timestep restriction in the simulation that scales like 
${\rm \Delta} x^2$ (Heath Reference Heath2002). We therefore choose a candidate timestep to satisfy 
${\rm \Delta} t = \zeta {\rm \Delta} x$ for a dimensionless constant 
$\zeta$, based on satisfying Courant–Friedrichs–Lewy conditions (Courant, Friedrichs & Lewy Reference Courant, Friedrichs and Lewy1967) for the advective terms in the liquid domain and gas layer.
The simulation outputs 
$N_f$ frames over the duration of the simulation. Thus, the time between frames is 
$t_f = t_{end}/N_f$. In general, an integer multiple of candidate timesteps will not exactly match 
$t_f$, so that 
$c{\rm \Delta} t < t_f < (c+1){\rm \Delta} t$. Because of this, the actual timestep is slightly adjusted to 
${\rm \Delta} t' = t_f/(c+1)$ and 
$c+1$ timesteps are taken between frames. Several examples demonstrating the performance of the code are provided in Appendix B.
4. Results and discussion
4.1. Initial dynamics and comparison with potential flow theory
To validate our approach, we first simulate the initial dynamics and we compare with the scaling results of potential flow theory that were introduced in § 2.4 for the height of the stagnation point 
$H^*$ (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012). We use the baseline physical parameters given in table 2. Since we only want to resolve the initial deceleration of the drop, and not the formation of the thin gas layer, we use the computational parameters in the third column of table 3, which feature a relatively coarse numerical grid of size 
$2048\times 256$. To compute 
$H^*$, we examine the sequence of height profiles from the simulation and find the time 
$t^*$ when 
$h_t(0,t^*)=0$, from which 
$H^* = h(0,t^*)$.
Baseline choices for the physical parameters used in the simulations. These parameters are used throughout the paper, with modifications to them noted in the text.
Baseline choices for the simulation parameters, which are all non-dimensional. The first set of parameters are for computing the initial dynamics and value of 
$H^*$ in § 4.1. The second two sets are for the lift-off calculations in all other sections. The procedure for connecting the non-dimensional variables 
$\tilde {L}$, 
$\tilde {H}_0$ and 
$\tilde {t}_{end}$ to physical values is described in § 2.4.
As described in § 2.4, at low initial velocities the flow in the gas layer can be treated as incompressible, and the scaling result 
$H^* = R \, St^{2/3}$ can be derived. Mani et al. (Reference Mani, Mandre and Brenner2010) extend this analysis to look at higher initial velocities, where compressibility of the gas becomes important. This happens at a height of
\begin{equation} H_* = R(\mu_g V/R P_0)^{1/2}, \end{equation}
where the subscript ‘
$*$’ is used to distinguish from the stagnation height. By modelling the subsequent drop deceleration below 
$H_*$ assuming gas compressibility, Mani et al. (Reference Mani, Mandre and Brenner2010) derive the result
\begin{equation} H^*_{comp} = R \, St^{2/3} \epsilon^{(2-\gamma)/(2\gamma-1)} \end{equation}for the stagnation height.
We first performed sequences of simulations using low initial velocities over the range from 
$V=0.15\ \text {m}\ \text {s}^{-1}$ to 
$V=1.35\ \text {m}\ \text {s}^{-1}$, using four different liquid viscosities from 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$ to 
$\nu _l = 300\ {\rm mm}^2\ {\rm s}^{-1}$. We also ran two sets of simulations with the surface tension set to half and zero its baseline value. Figure 3(a) shows a plot of 
$H^*/(R\, St^{2/3})$ as a function of 
$\epsilon ^{-1}$, with these six sets of data points in the left half of the domain where 
$\epsilon ^{-1}<1$ and gas compressibility is not important. We see that 
$H^*/(R\, St^{2/3})$ is approximately constant and is in agreement with (2.12). To examine the trends more clearly, figure 3(b) shows a zoomed-in plot of the same data. As expected, the best agreement is achieved for the smallest value of 
$\nu _l$ and for zero surface tension. For 
$\epsilon ^{-1}<1$, the values of 
$H^*/(R\, St^{2/3})$ are slightly lower for the cases of larger 
$\nu _l$ and 
$\sigma$. The decrease in height of the stagnation point as 
$\nu _l$ increases is qualitatively consistent with the experimental observations of Langley et al. (Reference Langley, Li and Thoroddsen2017), and hints at the limits of the potential flow theory (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010). By performing experiments with viscosities of up to 
$2 \times 10^6\ \text {mm}^2\ \text {s}^{-1}$, Langley et al. (Reference Langley, Li and Thoroddsen2017) are able to demonstrate that the centreline air thickness (closely related in definition to 
$H^*$) scales like 
$\mu _l^{-1/9}$. However, since we use inviscid boundary conditions on the top and sides of the fluid domain, we are unable to simulate large viscosities to compare to this result. Modifying our boundary conditions to work with large viscosities remains a subject for future work.
(a) Plot of rescaled initial stagnation height 
$H^*/(R\, St^{2/3})$ as a function of the inverse of the dimensionless compressibility parameter 
$\epsilon = P_0 R \, St^{4/3}/\mu _g V$, for a range of different liquid viscosities 
$\nu _l$, surface tensions 
$\sigma$ and gas parameters 
$\gamma$. Unless otherwise stated, baseline parameters from tables 2 and 3 are used. By default, the drop starts from a height 
$H_0$ scaled according to (2.13). For 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$, to account for compressibility effects, data from two simulation sequences that use a fixed initial height (FIH) based on substituting 
$V=2\ \text {m}\ \text {s}^{-1}$ into (2.13) are shown. (b) Zoomed-in plot of the same data, showing the region bounded by the dotted grey rectangle in panel (a).
We also performed two sequences of simulations with higher initial velocities from 
$V=1\ \text {m}\ \text {s}^{-1}$ to 
$V=18\ \text {m}\ \text {s}^{-1}$ using 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$, to examine the compressible regime. However, we note from (2.13) that our usual choice for initial height scales according to 
$V^{-2/3}$, whereas from (4.1), the height where compressibility is important scales according to 
$V^{1/2}$. To fully capture the compressible effects, the drop must start higher than 
$H_*$, and, thus, the usual initialization procedure is problematic for large 
$V$. For these simulations, we therefore set 
$H_0$ at a fixed value based on using 
$V=2\ \text {m}\ \text {s}^{-1}$ in (2.13). We ran two sequences of simulations using the usual choice of 
$\gamma =1.4$, and another with 
$\gamma =1$; as shown in figure 3, these results are in very good agreement with the scalings of 
$\epsilon ^{1/3}$ and 
$\epsilon$, respectively, that are expected from (4.2).
4.2. Evolution of dynamics and effect of liquid viscosity
With the initial dynamics validated, we now turn attention to simulating the continued spreading of the liquid drop on a layer of gas, and the effect of viscosity on the evolution of the liquid–gas interface. We use the same baseline physical parameters given in table 2, but since we must now accurately resolve the gas layer as it becomes very thin, we increase the simulation resolution as shown in the fourth and fifth columns of table 3. Furthermore, since the deviation of the interface from the solid surface happens later for high viscosities (Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014b), we use a larger domain and longer duration when 
$\nu _l > 20\ {\rm mm}^2\ {\rm s}^{-1}$, as indicated in the fifth column of table 3. We begin by using the baseline initial drop velocity of 
$V=0.45\ \text {m}\ \text {s}^{-1}$ to match the experiments of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b).
Figure 4 shows the height profiles for three different simulations using liquid viscosities of 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$, 
$\nu _l = 32\ {\rm mm}^2\ {\rm s}^{-1}$, and 
$\nu _l = 100\ {\rm mm}^2\ {\rm s}^{-1}$. Panels (a)–(c) show a large-scale view, where the initial parabolic profile approaches the surface, begins to decelerate and deforms to create the dimple. After the drop continues to fall, a thin layer of gas from 
$x\gtrapprox 200\ \mathrm {\mu }\text {m}$ is created and spreads outward. In all three simulations, the height profiles have a front that sweeps outward as more liquid approaches the surface.
Profiles of the height 
$h(x,t)$ of the gas layer at intervals spaced 
$6.004\ \mathrm {\mu }\text {s}$ apart, corresponding to an integer multiple of the frame output interval 
$t_f$ described in § 3.4, for liquid viscosities of (a) 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$, (b) 
$\nu _l=32\ {\rm mm}^2\ {\rm s}^{-1}$ and (c) 
$\nu _l=100\ {\rm mm}^2\ {\rm s}^{-1}$. All other simulation parameters follow the baseline values in tables 2 and 3. Panels (d)–( f) show the same data as (a)–(c), respectively, but with a smaller range of 
$h$ to highlight the lift-off behaviour. For each profile, the global minimum, which follows the leading tip, is also plotted on the curves; once the global minimum is no longer at the leading tip, it is no longer plotted. The dashed box in panel (d) marks a further zoomed-in region shown in figure 6.
Figure 5 shows snapshots of the pressure and vorticity in the liquid domain for the simulation with 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$. At 
$t=24.02\ \mathrm {\mu }\text {s}$, the drop is still approaching the surface. The pressure builds up close to 
$x=0$. Since the liquid near 
$x=0$ is decelerated faster, this creates a region of negative vorticity over 
$0 \lessapprox x \lessapprox 250\ \mathrm {\mu }\text {m}$. By 
$t=72.05\ \mathrm {\mu }\text {s}$, the thin layer of gas has formed, and the position of the front from figure 4(a) is marked with a small triangle. There is a region of large positive pressure behind the front, a small region of negative pressure ahead of it. The contour of zero vorticity follows the front as it moves outward to 
$t=96.07\ \mathrm {\mu }\text {s}$. The pressure fields are similar to those reported by Philippi et al. (Reference Philippi, Lagrée and Antkowiak2016), who also compute the pressure in the interior of the drop, using a different simulation framework (Popinet Reference Popinet2003, Reference Popinet2009; Lagrée et al. Reference Lagrée, Staron and Popinet2011).
Snapshots of pressure, 
$p$, (left) and vorticity, 
$\omega =[\boldsymbol {\nabla } \times \boldsymbol {u}]_3$, for a portion of the liquid domain for a simulation with liquid viscosity 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$, where all other parameters follow the baseline values in tables 2 and 3. The field values exhibit large variations in scale and also switch sign, so a nonlinear mapping 
$f(\alpha ) = ({1}/{\log 10}) \sinh ^{-1} ({\alpha }/{2})$ is used to create the colour maps. The thick black lines indicate the zero contour. The thin black lines show contours for 
$\pm 10^{n/2}\ \text {Pa}$ for pressure and 
$\pm 10^n\ \text{s}^{-1}$ for vorticity, where 
$n\in \mathbb {N}_0$. For 
$t=72.05\ \mathrm {\mu }\text {s}$ and 
$t=96.07\ \mathrm {\mu }\text {s}$, the position of the front (given by the local minimum of the profiles in figure 4a,d) is marked with a triangle.
Figure 4(d–f) shows zoomed-in plots of the height profiles in the thin layer for the three simulations. Behind the front, the height profiles align on top of each other, and trace out relatively stable envelopes, appearing as thick blue lines in figure 4(d–f). In all three cases the envelopes initially slope downward before curving upward, indicating the deviation of the liquid–gas interface away from the solid surface, prior to contact between the solid and the liquid. This is the lift-off phenomenon (Kolinski et al. Reference Kolinski, Mahadevan and Rubinstein2014b). Lift-off occurs more quickly for lower viscosities, and the envelope slopes upward faster. These results are in close agreement with the experimental results of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b).
We now quantify the lift-off time and position. While the envelopes in figure 4(d–f) are relatively well defined, the height profiles shift slightly in over time, making it difficult to precisely define the moment of lift-off. However, close inspection of figure 4(d–f) shows that in all cases the height profiles have a leading tip that dips slightly below the envelope that forms. We found this to be a well-defined feature that occurs universally across all of our simulations. The tip position can be determined as the global minimum of the height profile at each time, and provides a way to precisely define when lift-off occurs, at the moment when the leading tip starts to rise.
Figure 6(a) shows a further zoomed-in plot of the height profiles during lift-off for the case of 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$. Once the leading tip has risen sufficiently, then it is no longer the global minimum. However, since lift-off has always occurred by the time the leading tip is rising, this does not cause any difficulty in identifying the lift-off time. It is natural to consider whether the leading tip is a real feature or a numerical artifact that emerges from discretization error and limited resolution. To address this, figure 6(b) shows a further zoomed-in region, depicting the leading tip in detail. Here, the blue circles show individual simulation grid points. The grid spacing is substantially smaller than the width of the leading tip, indicating that it is a real feature. Further numerical tests of accuracy are provided in Appendix C.
Zoomed-in plots of the height of the gas layer at intervals spaced 
$1.2009\ \mathrm {\mu }\text {s}$ apart for liquid viscosity 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$, where all other parameters follow the baseline values in tables 2 and 3. Panel (a) shows the region marked by the dashed box in figure 4(d). For each profile, the global minimum, which follows the leading tip, is also plotted on the curves; once the global minimum is no longer at the leading tip, it is no longer plotted. Panel (b) shows the region marked by the dashed box in panel (a). In panel (b) the small blue circles indicate the computational grid. The grey dashed lines show the profiles with Gaussian smoothing applied, and the grey circles show the global minima of the smoothed lines.
The width of the tip is governed by surface tension. While many of our simulations use the baseline value of 
$\sigma = 0.072 \ \text {N}\ \text {m}^{-1}$ from table 2, we also consider reduced surface tension where the tip becomes sharper. In this case, there can be numerical difficulties in identifying the minimum due to per grid point variations in the height profile. To create a scheme for identifying the minimum across all simulations, we therefore smooth the height profile using a Gaussian kernel with length scale 
$1.2\ \mathrm {\mu }\text {m}$. This results in a minimal alteration in the leading tip position (figure 6b), but improves robustness when analysing the simulations. Numerically, the global minimum is found by identifying the local minima of the cubic interpolant of the smoothed profile, and selecting the smallest one.
Figures 7 and 8 show trajectories of the leading tip for a range of viscosities, from 
$\nu _l=2.5\ {\rm mm}^2\ {\rm s}^{-1}$ to 
$\nu _l=160\ {\rm mm}^2\ {\rm s}^{-1}$, for drop impact velocities of 
$V=0.45\ \text {m}\ \text {s}^{-1}$ and 
$V=0.9\ \text {m}\ \text {s}^{-1}$, respectively. The distance 
$x$ that characterises the horizontal extent of the dimple, as shown in figures 7(a) and 8(a), closely matches the experimental results of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) and Langley et al. (Reference Langley, Li and Thoroddsen2017), respectively. The height of the gas layer, 
$h \sim {O}(100)\ \text {nm}$, is also in a similar regime as the experimental observations of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) and Langley et al. (Reference Langley, Li and Thoroddsen2017), though larger in value than the 
$h$ observed by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b).
Trajectories of the global minimum of the drop profile in (a) the 
$(x,h)$ plane and (b) the 
$(t,h)$ plane for the baseline parameters with initial drop velocity 
$V=0.45\ \text {m}\ \text {s}^{-1}$, for a range of different liquid viscosities. All other parameters follow the baseline values in tables 2 and 3. For each trajectory, the filled circle indicates where the global minimum reaches its lowest point, which we define as when lift-off occurs. Each cross indicates when the global minimum no longer marks the leading tip. Trajectories for selected 
$\nu _l$ are labelled.
Trajectories of the global minimum of the drop profile in (a) the 
$(x,h)$ plane and (b) the 
$(t,h)$ plane for the baseline parameters with initial drop velocity 
$V=0.9\ \text {m}\ \text {s}^{-1}$, for a range of different liquid viscosities. All other parameters follow the baseline values in tables 2 and 3. For each trajectory, the filled circle indicates where the global minimum reaches its lowest point, which we define as when lift-off occurs. Each cross indicates when the global minimum no longer marks the leading tip. Trajectories for selected 
$\nu _l$ are labelled.
Two regimes are visible in figure 7. For 
$\nu _l \lessapprox 20\ {\rm mm}^2\ {\rm s}^{-1}$, increasing viscosity results in the trajectory reaching a lower value of 
$h$, and the lift-off position moving slightly outward. For 
$\nu _l \gtrapprox 20\ {\rm mm}^2\ {\rm s}^{-1}$, increasing viscosity results in the trajectory reaching a higher value of 
$h$, and the lift-off position moves outward substantially. Small undulations in the trajectories are visible for 
$\nu _l \gtrapprox 20\ {\rm mm}^2\ {\rm s}^{-1}$, which arise due to capillary waves in the height profile. This can have a substantial effect on the measured lift-off time, depending on which undulation achieves the global minimum. For example, there is a large difference between 
$\nu _l = 20\ {\rm mm}^2\ {\rm s}^{-1}$ and 
$\nu _l = 25\ {\rm mm}^2\ {\rm s}^{-1}$. In figure 8, for 
$V=0.9\ \text {m}\ \text {s}^{-1}$, the scale of 
$h$ is smaller. The observed reduction in film thickness with increasing 
$V$ is consistent with the experimental observations of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b). Considering the case of liquid–substrate contact being mediated by a reduced thickness of the air layer, our observations are also consistent with the transition from gliding over an air layer to contact with increasing 
$V$, as reported by Langley et al. (Reference Langley, Li and Thoroddsen2017). Figure 8 also shows a larger relative amplitude of the capillary waves as compared with the smaller impact velocity of figure 7. As a result, the capillary waves have a noticeable effect on lift-off times at lower viscosities, with a large difference in lift-off time between 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$ and 
$\nu _l =13\ {\rm mm}^2\ {\rm s}^{-1}$.
We now compare to the experimental finding of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) that the lift-off time is proportional to 
$\nu _l^{1/2}$. The lift-off time 
$\tau$ in this previous work is defined relative to a time origin of when the drop would reach 
$h=0$ in the absence of the surface. Here, since the initial height 
$h_0$ and velocity 
$V$ are known, we compute the time origin as 
$t_0=h_0/V$. By contrast, Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) were not able to directly view 
$h_0$, since the drop can only be observed once it enters an evanescent field of height 
$h_{ev}=1\ \mathrm {\mu }\text {m}$. Instead, they measured the first position and time 
$(h',r',t')$ when the drop enters the evanescent field, and assume that the drop is undeformed at this position, so that 
$h' = h_0 - Vt' + r'^2/(2R^2)$. Hence, 
$h_0$ can be calculated and, therefore, the time origin is 
$t^{exp}_0=h_0/V - t' + r'^2/(2R^2V)$.
Figure 9 shows the lift-off times 
$\tau$ as a function of viscosity for seven different values of initial velocity 
$V$. The data points from Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) are also plotted. Even though the experiments of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) are performed in a three-dimensional axisymmetric configuration, there is good quantitative agreement with the two-dimensional simulation results.
Lift-off time 
$\tau$ as a function of liquid viscosity 
$\nu _l$, for a range of initial drop velocities 
$V$. All other parameters follow the baseline values in tables 2 and 3. The experimental data are taken from figure 4(a) of the paper by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b), which used 
$V=0.45\ \text {m}\ \text {s}^{-1}$. Three data points are omitted from the plot for low 
$V$ and low 
$\nu _l$ due to physical difficulties with running the simulation; see § 4.3.
Overall, the results are consistent with the 
$\nu _l^{1/2}$ scaling result, but the additional precision provided by the simulations reveals a more complicated relationship between 
$\nu _l$ and 
$\tau$. For each value of 
$V$, the data points exhibit some fluctuations, which are due to the undulations visible in figure 7 where the global minimum defining the lift-off time may abruptly change. While a 
$\nu _l^{1/2}$ scaling appears consistent with the high impact velocities where 
$V\gtrapprox 0.75\ \text {m}\ \text {s}^{-1}$, the data for low impact velocities looks a better fit to 
$\nu _l^\eta$, where there are two different values of 
$\eta$ with a change at 
$\nu _l \approx 20\ {\rm mm}^2\ {\rm s}^{-1}$. This is also consistent with the two different types of behaviour observed in figure 7.
The slope of data points for small values of 
$\nu _l$ in figure 9 is strongly affected by the choice of time origin, and may affect the conclusions about the relevant exponents. To investigate this further, we replotted the data using 
$t_0^{exp}$ as the time origin; this resulted in a small shift upwards of the data (since 
$t_0^{exp}< t_0$ in all cases), but did not affect the overall patterns. As a third approach, we defined a time origin 
$t_0^{dat}$ directly from the data, by finding the best fit to the model 
$\tau _{dat} = t-t_0^{dat} = \gamma (\nu _l/\nu _{sc})^\alpha$ for the three free parameters 
$(t_0^{dat}, \gamma, \alpha )$. Here 
$\nu _{sc}=20\ {\rm mm}^2\ {\rm s}^{-1}$ is chosen as an arbitrary viscosity scale. Specifically, for the data points 
$(\nu _{l,k},t_k)$, we minimized the residual
\begin{equation} S(t_0^{dat},\gamma,\alpha) =\frac{1}{2} \sum_k \left( \log \gamma + \alpha \log \frac{\nu_{l,k}}{\nu_{sc}} - \log(t_k-t_0^{dat})\right)^2. \end{equation}
We used the L-BFGS-B algorithm for bound-constrained nonlinear optimization (Byrd et al. Reference Byrd, Lu, Nocedal and Zhu1995), and enforced 
$t_0^{dat}<0.999 t_{min}$ and 
$\alpha >0$, where 
$t_{min} = \min _k \{t_k\}$. We started the minimization using multiple initial guesses with 
$t_0^{dat} \in [0.7t_{min},0.9t_{min}]$ and 
$\alpha \in [0.4,1.2]$. For each pair 
$(t_0^{dat},\alpha )$, the value of 
$\gamma$ is set so that the bracketed expression in (4.3) is zero for the 
$\nu _l=20\ {\rm mm}^2\ {\rm s}^{-1}$ data point, so that the initial guess should be close to the minimum of 
$S$.
Figure 10 shows a replotting of the data relative to this time origin. The data for 
$V=0.15\ \text {m}\ \text {s}^{-1}$ is omitted from this plot since the lift-off times are non-monotonic for the smallest values of 
$\nu _l$. For each other value of 
$V$, the minimization procedure converges to a single unique solution for all initial guesses. With this definition of time origin, all of the data are more consistent with a linear scaling relationship as opposed to the square-root relationship in figure 9. Panels (b)–(d) show the values of the fitted parameters, and demonstrate that 
$\alpha \in [1.02,1.28]$ in all cases. The average residual over the six different values of 
$V$ is 
$\bar {S}=0.2536$. If the exponent is constrained to 
$\alpha =\frac {1}{2}$ to match the exponent of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b) then the average residual increases to 
$\bar {S}=0.3855$ and the data points exhibit an upward curve in a log–log plot that systematically deviates from a power law. Constraining the exponent to 
$\alpha =1$ results in 
$\bar {S}=0.2805$ and a better fit to the model that is similar to the three-parameter fit shown in figure 10. See the supplementary information for additional discussion and parameter fits. These results highlight that any theoretical analysis of the relationship between 
$\tau$ and 
$\nu _l$ would need to take into account the sensitivity of defining the time origin.
(a) Lift-off time 
$\tau _{dat}$ as a function of liquid viscosity 
$\nu _l$, for a range of initial drop velocities 
$V$, using the alternative time origin definition based on by minimizing the three-parameter residual function in (4.3). All other parameters follow the baseline values in tables 2 and 3. (b–d) Best fits of the parameters 
$t_0^{dat}- t_0$, 
$\gamma$, and 
$\alpha$ for different initial drop velocities.
Figure 11(a) shows lift-off times (relative to 
$t_0$) for three different drop radii: the original value of 
$R=1.5\ \text {mm}$, and as well as half and double this value. Increasing the radius increases the lift-off time, similar to the effect of lowering velocity in figure 9.
Lift-off time as a function of liquid viscosity, for a range of (a) different drop radii and (b) different surface tension values. All other parameters follow the baseline values in tables 2 and 3. The data point for 
$(\nu _l,R)=(2.5\ {\rm mm}^2\ {\rm s}^{-1}, 0.75\ \text {mm})$ is omitted due to physical difficulties with running the simulation (§ 4.3). The data points for 
$\sigma =0.0072\ \text {N}\ \text {m}^{-1}$ and 
$\nu _l\ge 100\ {\rm mm}^2\ {\rm s}^{-1}$ are omitted because lift-off does not occur over the simulation duration.
4.3. The role of surface tension
The simulations allow us to change the surface tension in ways that would be difficult to do experimentally, to investigate the importance of this physical effect. Changing the surface tension allows us to suppress or accentuate the capillary waves in the liquid–gas interface, as shown in figure 7, to investigate the effect of surface tension on the phenomenon of lift-off. Following the same procedure as in § 4.2, we calculated the lift-off times for 
$\sigma =0.0072\ \text {N}\ \text {m}^{-1}$, 
$\sigma = 0.036\ \text {N}\ \text {m}^{-1}$ and the 
$\sigma = 0.144\ \text {N}\ \text {m}^{-1}$, corresponding to a tenth, half and double the usual surface tension values, respectively. Figure 11(b) shows that as surface tension is reduced, the lift-off times increase markedly. For the case of 
$\sigma =0.0072\ \text {N}\ \text {m}^{-1}$, lift-off is eliminated for large values of 
$\nu _l$, with the leading tip continuing to decrease over the course of the simulation.
Figure 11(b) strongly suggests that surface tension is important in creating lift-off. Reducing from 
$\sigma = 0.036\ \text {N}\ \text {m}^{-1}$ to 
$\sigma =0.0072\ \text {N}\ \text {m}^{-1}$ almost doubles the lift-off times in most cases, and one can ask whether lift-off will be completely eliminated in the limit as surface tension vanishes. To examine this, we ran a sequence of simulations with 
$\sigma =0$, with several representative examples for 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$, 
$\nu _l = 32\ {\rm mm}^2\ {\rm s}^{-1}$ and 
$\nu _l=100\ {\rm mm}^2\ {\rm s}^{-1}$ shown in figure 12. This is a difficult limit to probe in our simulations, since as discussed for figure 6, surface tension regularizes the leading tip. Figure 13 shows close-ups of the profiles for 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$, indicating that the leading tip becomes as sharp as a single grid point, and the minima of the profiles can fluctuate non-monotonically depending on exactly how the tip aligns with the computational grid. However, in this case the profile smoothing procedure introduced in § 4.2 is sufficient to extract smooth leading tip trajectories.
Profiles of the height of the gas layer at intervals spaced 
$6.004\ \mathrm {\mu }\text {s}$ apart for liquid viscosities of (a) 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$, (b) 
$\nu _l=32\ {\rm mm}^2\ {\rm s}^{-1}$ and (c) 
$\nu _l=100\ {\rm mm}^2\ {\rm s}^{-1}$ using zero surface tension. All other parameters follow the baseline values in tables 2 and 3. Panels (d)–( f) show the same data as (a)–(c), respectively, but with a smaller range of 
$h$ to highlight that no lift-off occurs in this case. For each profile, the global minimum, which follows the leading tip, is also plotted on the curves. The dashed box in panel (d) marks a further zoomed-in region shown in figure 13.
Zoomed-in plots of the height of the gas layer at intervals spaced 
$1.2009\ \mathrm {\mu }\text {s}$ apart for liquid viscosities of 
$\nu _l = 10\ {\rm mm}^2\ {\rm s}^{-1}$ using zero surface tension. All other parameters follow the baseline values in tables 2 and 3. Panel (a) shows the region marked by the dashed box in figure 12(d). For each profile, the global minimum, which follows the leading tip, is also plotted on the curves. Panel (b) shows the region marked by the dashed box in panel (a). In panel (b) the small blue circles indicate the computational grid. The grey dashed lines show the profiles with Gaussian smoothing applied, and the grey circles show the global minima of the smoothed lines.
The very sharp leading tip causes another difficulty in the simulations: as the tip is advected across the computational grid, there will be an effective numerical diffusion, which will act as though a small surface tension has been imposed. This is an important issue since figure 11(b) already confirms that small surface tensions can considerably alter the behaviour. To test this, we compared simulations with the baseline parameters to those on finer grids (Appendix C). Unlike the case for finite surface tension, the simulations on finer grids are noticeably different, with the profiles reaching lower heights and the leading tip not curving up as rapidly. Since liquid viscosity also regularizes the evolution of the height profile, these discrepancies are more significant when 
$\nu _l$ is small.
We calculated the trajectories of the leading tip for a range of values for liquid viscosity, 
$\nu _l$. For 
$13\ {\rm mm}^2\ {\rm s}^{-1} <\nu _l \le 40\ {\rm mm}^2\ {\rm s}^{-1}$, we switched to a larger computational grid of size 
$8192\times 1536$ and, for 
$\nu _l \le 13\ {\rm mm}^2\ {\rm s}^{-1}$, we switched to a very large grid of size 
$12\ 288 \times 2304$. For 
$\nu _l \le 44\ {\rm mm}^2\ {\rm s}^{-1}$, we were not able to simulate on a large enough grid to adequately resolve the numerical diffusion and, hence, results in this range are omitted. Figure 14 shows the trajectories in both the 
$(x,h)$ and 
$(t,h)$ planes. In all cases, the leading tips continue to decrease. While our results only examine one particular set of physical parameters (from table 2) our results strongly suggest that surface tension is required for lift-off to occur. Moreover, we observe that the thickness of the gas layer increases with 
$\nu _l$. This is consistent with the experimental findings of Langley et al. (Reference Langley, Li and Thoroddsen2017), who observe a transition from drop–substrate contact to gliding of the drop over a thin air film as 
$\nu _l$ increases.
Trajectories of the global minimum of the drop profile in (a) the 
$(x,h)$ plane and (b) the 
$(t,h)$ plane for a range of different liquid viscosities, using the baseline parameters in tables 2 and 3 and zero surface tension. Simulations with 
$13\ {\rm mm}^2\ {\rm s}^{-1} <\nu _l \le 40\ {\rm mm}^2\ {\rm s}^{-1}$ use a grid of size 
$8192\times 1536$, and simulations with 
$\nu _l \le 13\ {\rm mm}^2\ {\rm s}^{-1}$ use a grid of size 
$12\ 288 \times 2304$. Trajectories for selected 
$\nu _l$ are labelled.
We also found a regime where the effect of surface tension can qualitatively affect the lift-off behaviour. Figure 15(a) shows the height profiles for a simulation with the baseline value of surface tension of 
$\sigma =0.072\ \text {N}\ \text {m}^{-1}$ in the regime of low initial velocity, 
$V=0.3\ \text {m}\ \text {s}^{-1}$ and low viscosity, 
$\nu _l = 6.5\ {\rm mm}^2\ {\rm s}^{-1}$. In this regime, prominent capillary waves are generated outside the thin gas layer, which grow larger as time progresses. Figure 15(b) shows a zoomed-in plot of the profiles in the thin gas layer. Lift-off appears to occur at 
$x \approx 280 \ \mathrm {\mu }\text {m}$ and the leading tip starts to rise. However, the influence of the capillary wave causes the tip to move downward again, ultimately dipping below the previous minimum height. It is possible that the tip may move upward again at a later point, but we were unable to track the behaviour further. The sharp features visible in figure 15(a) (e.g. at 
$(x,h)\approx (700\ \mathrm {\mu }\text {m},100\ \mathrm {\mu }\text {m})$) cause the simulation to terminate early, since the Newton–Raphson iterations can no longer be solved to an acceptable level of accuracy.
(a) Profiles of the height of the gas layer at intervals spaced 
$11.80\ \mathrm {\mu }\text {s}$ apart for a liquid viscosity of 
$\nu _l=6.5\ {\rm mm}^2\ {\rm s}^{-1}$ and initial drop velocity of 
$V=0.3\ \text {m}\ \text {s}^{-1}$, showing the development of capillary waves. All other parameters use the baseline values in tables 2 and 3. (b) Zoomed-in plot of the same data. The global minima of the curves are plotted when they indicate the leading tip.
4.4. Effect of compressibility in the gas
We performed simulations to examine the lift-off behaviour in the regime of high initial drop velocity when compressibility in the gas becomes important. We chose 
$V=2\ \text {m}\ \text {s}^{-1}$, which is substantially into the compressible regime of figure 3. We focused on 
$\nu _l=32\ {\rm mm}^2\ {\rm s}^{-1}$, and due to the faster dynamics in the regime, we halved the timestep multiplier from its baseline value to 
$\zeta = 4\times 10^{-3}$. Figure 16(a) shows that the behaviour is qualitatively different in this case, with the central dimple noticeably rebounding while the drop profile spreads outward. Figure 16(b) shows a zoomed-in region of the thin gas layer. Other than 
$V$, all physical simulation parameters are identical to those used in figure 4(e), yet the results are considerably different with a much thinner gas layer, and visually, they are a closer match to the simulations without surface tension. For these simulations, we were therefore not able to identify the lift-off time. Figure 16(c) shows plots of the centreline height 
$h(0,t)$ for four different values of 
$V$ indicating that the relative rebound of the central dimple is much larger for higher drop impact velocities.
(a) Profiles of the height of the gas layer spaced 
$0.4998\ \mathrm {\mu }\text {s}$ apart for a liquid viscosity of 
$\nu _l=32\ {\rm mm}^2\ {\rm s}^{-1}$ and an initial drop velocity of 
$V=2\ \text {m}\ \text {s}^{-1}$, where the gas compressibility becomes important. The timestep multiplier is set to 
$\zeta = 4\times 10^{-3}$. All other parameters use the baseline values in tables 2 and 3. (b) Zoomed-in plot of the same data. (c) Semi-log plot of the centreline height 
$h(0,t)$ for four different drop velocities when 
$\nu _l=32\ {\rm mm}^2\ {\rm s}^{-1}$, with all other parameters using the baseline values.
5. Conclusion
In this paper we investigated the viscous effects in the early stages of drop impact on a surface. We coupled a Navier–Stokes solver to model flow in the interior of the drop with a partial differential equation to model the pressure and height in the thin gas layer between the drop and the surface. We demonstrated that our simulations are consistent with previous work using potential flow theory where flow in the liquid is assumed incompressible (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012). However, our simulations allow us to go beyond this previous work and investigate viscous effects. We showed that at low initial drop velocities, viscosity plays a weak role in the deceleration of the drop, and the height 
$H^*$ at which it reaches a stagnation point.
Using our simulations, we are able to recreate the lift-off phenomenon that was experimentally reported by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b). We have therefore demonstrated that with the reduced model described in § 2, our simulations are able to capture lift-off. Our numerical results for the lift-off time 
$\tau$ are consistent with the 
$\nu _l^{1/2}$ scaling relationship found by Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b). However, the additional precision afforded by the simulations indicates that the precise relationship between 
$\tau$ and 
$\nu _l$ is more complex, due to the effects of capillary waves and the sensitivity of the results to the definition of the time origin. The simulation allows us to probe conditions that would be difficult to observe experimentally. Using this capability, our results provide strong evidence that surface tension is necessary for lift-off to occur. This is consistent with the numerical study of Duchemin & Josserand (Reference Duchemin and Josserand2011), who found that surface tension was required for the formation of a thin jet skating above the surface. In their study, the absence of surface tension led to a finite-time singularity where the liquid–gas interface touches the substrate, whereas in our case the interface asympotically approaches the substrate. However, there are several differences between their simulations and ours, such as the usage of axisymmetric coordinates, and a different way to represent the liquid–gas interface and solve for the relevant physical variables. It is therefore difficult to make an exact comparison.
There are a number of possible next steps. The simulations provide a detailed view of the lift-off phenomenon, and all aspects (e.g. stress, velocity, vorticity) can be calculated. The results may therefore guide theoretical analyses of lift-off, and may allow scaling relationships similar to those presented by Mandre et al. (Reference Mandre, Mani and Brenner2009), Mandre & Brenner (Reference Mandre and Brenner2012) and Mani et al. (Reference Mani, Mandre and Brenner2010) to be derived. As part of our study, we have released a complete high-performance open source code that can examine lift-off across a wide range of configurations.
We opted to use a fixed-grid simulation for simplicity in coupling the liquid domain and gas layer together, and as described in Appendices B and C we are able to obtain accurate simulation results in a reasonable timeframe. However, it is likely that the behaviour at the bottom of the liquid domain, close to the liquid–gas interface and near the leading tip, dominates. Because of this, the simulation is a good candidate for use with adaptive mesh refinement (Berger & Oliger Reference Berger and Oliger1984; Guittet, Theillard & Gibou Reference Guittet, Theillard and Gibou2015), where the liquid–gas interface and the region around the leading tip could use a finer mesh. Adaptive mesh refinement has already been used for full-scale simulations of droplet impact, such as those using the Gerris/Basilisk software (Popinet Reference Popinet2003, Reference Popinet2009; Lagrée et al. Reference Lagrée, Staron and Popinet2011; Popinet Reference Popinet2015; Philippi et al. Reference Philippi, Lagrée and Antkowiak2016), but it could also be a useful avenue to explore for our reduced model. It would result in substantial computational savings, although the liquid–gas coupling would become more complicated and numerical errors would be harder to quantify.
The simulation could also be generalized to use cylindrical axisymmetric coordinates. The overall numerical approach would stay the same, but additional radial factors would have to be incorporated throughout the simulation. The Navier–Stokes solver that we employ has already been demonstrated to work in axisymmetric coordinates (Yu et al. Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007), but the routines that involve the gas layer would require modifications. For example, the kernel in (2.10) would need to be modified. While the current simulation is already in good agreement with the experimental results of Kolinski et al. (Reference Kolinski, Mahadevan and Rubinstein2014b), an axisymmetric solver would allow for a near-perfect comparison. This would, for example, help elucidate further the precise relationship between lift-off time 
$\tau$ and liquid viscosity 
$\nu _l$.
Experiments by Thoroddsen and coworkers have used interferometry to examine the evolution of the gas layer across a full two-dimensional surface (Langley et al. Reference Langley, Li and Thoroddsen2017). Their results show a number of interesting effects beyond the scope of our current model, such as a breakage of rotational symmetry and the formation of ruptures in the air film (Langley et al. Reference Langley, Li and Thoroddsen2017; Li et al. Reference Li, Langley, Tian, Hicks and Thoroddsen2017). They have also examined the case of nano-rough surfaces (Langley et al. Reference Langley, Li, Vakarelski and Thoroddsen2018). To connect with this work, our model could be generalized to a full three-dimensional simulation of the liquid and two-dimensional simulation of the gas layer. This would be substantially more computationally challenging, and would be a good candidate for using parallel computing and adaptive mesh refinement (Zhang et al. Reference Zhang2019).
Funding
C.H.R. was partially supported by the National Science Foundation under grant no. DMS-1753203. C.H.R. was partially supported by the Applied Mathematics Program of the U.S. DOE Office of Science Advanced Scientific Computing Research under contract no. DE-AC02-05CH11231.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The C++ code to generate all of these results is provided in the GitHub repository at https://github.com/chr1shr/vdropimpact. The code makes use of the TGMG library, which is provided in a separate GitHub repository at https://github.com/chr1shr/tgmg.
Appendix A. Calculations for the governing equation in the gas layer
A.1. Derivation of the gas layer pressure update equation
Substituting for 
$\rho _g$ from the equation of state, (2.5), into the lubrication equation for the pressure in the gas layer, (2.4), we get
\begin{equation} 12 \mu_g \left(\left(\frac{\rho_0}{p_0^{1/\gamma}} p_g^{1/\gamma} \right)h\right)_t = \left(\left(\frac{\rho_0}{p_0^{1/\gamma}} p_g^{1/\gamma} \right) h^3 p_{g,x}\right)_x. \end{equation}
Dividing both sides by 
${\rho _0}/{p_0^{1/\gamma }}$,
\begin{equation} 12 \mu_g ( p_g^{1/\gamma} h)_t = (p_g^{1/\gamma} h^3 p_{g,x})_x. \end{equation}Using the product rule to expand the brackets,
\begin{equation} 12 \mu_g \left( \frac{1}{\gamma}p_g^{1/\gamma-1} p_{g,t}h + p_g^{1/\gamma} h_t \right) = \frac{1}{\gamma}p_g^{1/\gamma-1}p_{g,x} h^3 p_{g,x} +p_g^{1/\gamma} 3 h^2 h_x p_{g,x} +p_g^{1/\gamma} h^3 p_{g,xx}. \end{equation}
Dividing both sides by 
$p_g^{1/\gamma -1}$,
\begin{equation} 12 \mu \left( \frac{1}{\gamma} p_{g,t}h + p_g h_t \right) = \frac{1}{\gamma} p_{g,x} h^3 p_{g,x} +p_g 3 h^2 h_x p_{g,x} +p_g h^3 p_{g,xx}. \end{equation}
Dropping the subscript 
$g$ for gas yields
\begin{equation} 12 \mu \left( \frac{1}{\gamma} p_{t}h + p h_t \right) = \frac{1}{\gamma} p_{x} h^3 p_{x} +p 3 h^2 h_x p_{x} +p h^3 p_{xx}, \end{equation}which can be rearranged to give (3.4).
A.2. Numerical solution of the gas layer pressure
We now describe how to solve (3.5) in the main text for updating the gas layer pressure using the Newton–Raphson method. Let 
$P^k$ be a vector containing the pressure estimates at the 
$k$th Newton–Raphson iteration. As an initial estimate, we set 
$P^0$ to be the pressure values from the 
$n$th timestep. We then write (3.5) as a nonlinear system 
$F(P)=0$, where the 
$i$th component of 
$F$ is given by the left-hand side minus the right-hand side of (3.5) evaluated at the 
$i$th grid point. An improved estimate 
$P^{k+1}$ for the pressure is given in terms of 
$P^k$ by
\begin{equation} J_F(P^k) (P^{k+1} - P^k) ={-}F(P^k), \end{equation}
where 
$J_F(P^k)$ is the Jacobian of 
$F$, which has components
\begin{equation} J^{k+1}_{F,ij} = \frac{\partial F_i}{\partial P^{k+1}_j}.\end{equation}
Since the finite-difference stencils in (3.5) only involve adjacent grid points, 
$J_F$ is a tridiagonal system. It can be written as
\begin{equation} J_F = \left( \begin{array}{ccccc} b_0 & c_0 & & & \\ a_1 & b_1 & c_1 & & \\ & a_2 & b_2 & c_2 & \\ & & \ddots & \ddots & \ddots \end{array} \right), \end{equation}where the terms are given by
\begin{align} a_i &= \frac{\partial F}{\partial p_{i-1}^{n+1}} = \frac{1}{\gamma} \frac{\bar{h}^3}{{\rm \Delta} x} \left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2{\rm \Delta} x}\right) - \frac{3 \bar{h}^2 \bar{h}_x }{2 {\rm \Delta} x} p^{n+1}_i + \frac{\bar{h}^3}{{\rm \Delta} x^2}{p^{n+1}_i}, \end{align}
\begin{align} b_i &= \frac{\partial F}{\partial p_{i}^{n+1}} ={-}12\frac{\mu}{\gamma}\frac{\bar{h}}{{\rm \Delta} t} (1) - 12 \mu \bar{h}_t -\left[\frac{3 \bar{h}^2 \bar{h}_x }{2 {\rm \Delta} x} \left({p^{n+1}_{i+1}-p^{n+1}_{i-1}}\right)\right. \nonumber\\ &= \left.+ \frac{\bar{h}^3}{{\rm \Delta} x^2} \left(p_{i+1}^{n+1} - 2 p_{i}^{n+1} + p_{i-1}^{n+1}\right) + \frac{\bar{h}^3}{{\rm \Delta} x^2} p_{i}^{n+1}({-}2) \right], \end{align}
\begin{align} c_i &= \frac{\partial F}{\partial p_{i+1}^{n+1}} ={-}\frac{1}{\gamma} \frac{\bar{h}^3}{{\rm \Delta} x} \left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2 {\rm \Delta} x}\right) - \frac{3 \bar{h}^2 \bar{h}_x }{2 {\rm \Delta} x} p^{n+1}_i - \frac{\bar{h}^3}{{\rm \Delta} x^2}{p^{n+1}_i}. \end{align}Solving (A6) can be done efficiently using LAPACK's tridiagonal solver dgtsv (Anderson et al. Reference Anderson1999). In most cases, since the pressure does not change by a large amount per timestep, fewer than five iterations are required in order to achieve numerical convergence.
Appendix B. Tests of the simulation performance
Table 4 contains statistics about the performance of the code for several simulations that were referenced in the main text. All tests were run using ten threads on an Ubuntu Linux computer with a ten-core 2.8 GHz Intel Core i9-10900 CPU, and the code was compiled with GCC version 10.3.
Performance statistics for several different simulations, compiled using GCC 10.3 on an Ubuntu Linux computer with an 2.8 GHz Intel Core i9-10900 CPU. Ten threads were used, and all simulations use the baseline parameter choices in tables 2 and 3 unless otherwise noted. The wall clock (WC) time is reported for each test, and the fraction of time on major components, such as the computation of boundary conditions (BCs), solving the MAC linear system and solving the FEM linear system, are reported.
The initial dynamics simulations only need to capture the large-scale deformation of the drop, and can therefore be run on relatively coarse computational grids. Consequently, a typical simulation takes approximately 40 min of wall clock time to run. A large portion of the computation time is spent on solving the linear systems with the multigrid method. This should be expected since the multigrid V-cycles involve repeated scans over the entire grid. Collectively, the four multigrid solves take up 44 % of the total computation time. The MAC and FEM linear systems take slightly more time and V-cycles, since apart from where Dirichlet boundary conditions are applied, these linear systems are only weakly diagonally dominant. By contrast, the linear systems for the implicit viscous term are strictly diagonally dominant. We note that the linear systems for the 
$u$ and 
$v$ velocity components are slightly different, since at 
$x=0$ we apply different boundary conditions, 
$(u,v_x)=(0,0)$. Hence, it is not possible to vectorize this system, and there is no substantial computational advantage over solving for the updates to 
$u$ and 
$v$ separately.
Calculating the boundary conditions according to (2.10) requires applying Simpson's rule along the 
$M$ grid points on the bottom edge. This must be done for the top edge of length 
$M$, and the right edge of length 
$N$, requiring 
${O}(M(M+N))$ work in total. This is a sizable amount of work and takes 3.4 % of the total computation time. Even though the gas layer involves several Newton steps and tridiagonal matrix solves, it only needs 
${O}(M)$ work and, therefore, takes up a minimal amount of the total computation time.
Table 4 also contains performance information for two lift-off simulations, which are run on larger computational grids to obtain high accuracy in the height of the gas layer. The wall clock time per timestep increases from 81 ms to 1.041 s. Based on a linear scaling with grid points, we would expect 81 ms to increase to 
$(81\ \text {ms})\frac {5120\times 960}{2048\times 256} = 760\ \text {ms}$. Thus, the performance is comparable, but slightly worse than, linear scaling, which may be due to reduced cache efficiency for a larger grid. Overall, the percentages spent on the different parts of the simulation are comparable, although the fraction spent on multigrid V-cycles increases slightly, and the fraction spent on the gas layer (which scales like 
${O}(M)$) decreases. In particular, more V-cycles are required to handle the implicit viscosity linear system when 
$\nu _l$ is larger.
Appendix C. Tests of the simulation accuracy
The simulation domain size, which is set via the non-dimensional parameter 
$\tilde {L}$ in table 3 and the domain aspect ratio 
$\beta$, may have an effect on the results. Since the boundary conditions on the top and right boundaries are based on an inviscid assumption, the domain size affects the extent to which viscosity is resolved. In addition, when solving for the pressure in the gas layer, the boundary condition of 
$p=P_0$ is imposed at 
$x=L$, and thus extending the domain affects the influence of this boundary condition on the simulation.
To test the sensitivity of the simulation results to the domain size, we performed simulations where either the horizontal dimension, vertical dimension or both dimensions were extended by a factor of 1.5. In each of these simulations the grid spacings 
${\rm \Delta} x = {\rm \Delta} y$ were kept the same, so that the horizontal extension increases the grid points from 5120 to 7680, and the vertical extension increases the grid points from 960 to 1440. Figure 17 shows the height profiles in the thin gas layer for the four simulations, for (a) 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$ and (b) 
$\nu _l=100\ {\rm mm}^2\ {\rm s}^{-1}$. The vertical extensions give visually indistinguishable curves, indicating that the vertical dimension is sufficiently large to resolve the viscous effects even for the larger value of 
$\nu _l$. The horizontal extensions create minor vertical shifts in the curves, suggesting that the pressure boundary condition has an effect on the results. However, these shifts are still acceptably small, and result in no substantial change to the position and time at which lift-off occurs.
Plots of the height profiles at 
$t=102.07\ \mathrm {\mu }\text {s}$ for simulations with (a) liquid viscosity 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$ and (b) liquid viscosity 
$\nu _l=100\ {\rm mm}^2\ {\rm s}^{-1}$. Baseline parameters from tables 2 and 3 are used although the original domain size for both values of 
$\nu _l$ uses 
$\tilde {L}=30$ and 
$\beta =\frac {16}{3}$. Results are also shown where the simulation domain is extended by a factor of 1.5 in either or both dimensions. In the extended simulations, the number of grid points is increased to keep the grid spacings 
${\rm \Delta} x={\rm \Delta} y$ the same.
We also examined the sensitivity of the results to the numerical grid size. Figure 18(a) shows height profiles for a simulation using the baseline parameters with liquid viscosity 
$\nu _l=10\ {\rm mm}^2\ {\rm s}^{-1}$ and surface tension 
$\sigma =0.072\ \text {N}\ \text {m}^{-1}$. Simulations with higher resolutions of 
$8192\times 1536$ and 
$12288\times 2304$ are also plotted, resulting in visually indistinguishable results. Figure 18(b) shows height profiles when the surface tension is set to zero and all other parameters are kept the same. In this case, as discussed in § 4.3, the leading tip becomes very sharp since the regularizing effect of surface tension is removed. As the tip moves across the grid, there will be a resolution-dependent numerical diffusion that will act as an effective small surface tension. Thus, in this case there is a small shift in the height profiles. While this remains small, it makes the precise behaviour of the leading tip difficult to resolve for small values of 
$\nu _l$, requiring larger grid sizes as presented in § 4.3.
Plots of the height profiles spaced 
$6.004\ \mathrm {\mu }\text {s}$ apart for three different resolutions using the baseline parameters from tables 2 and 3, and liquid viscosity 
$\nu _l =10\ {\rm mm}^2\ {\rm s}^{-1}$, with (a) surface tension 
$\sigma =0.072\ \text {N}\ \text {m}^{-1}$ and (b) zero surface tension.
Open access
