# Computing the viscous effect in early-time drop impact dynamics

Shruti Mishra1  Shmuel M. Rubinstein1,2  Chris H. Rycroft1,3
<EMAIL_ADDRESS>1John A. Paulson School of Engineering and Applied
Sciences, Harvard University, Cambridge, MA 02138, USA 2The Racah Institute of
Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
3Computational Research Division, Lawrence Berkeley Laboratory, Berkeley, CA
94720, USA

###### Abstract

The impact of a liquid drop on a solid surface involves many intertwined
physical effects, and is influenced by drop velocity, surface tension, ambient
pressure and liquid viscosity, among others. Experiments by Kolinski et al.
(2014b) show that the liquid–air interface begins to deviate away from the
solid surface even before contact. They found that the lift-off of the
interface starts at a critical time that scales with the square root of the
kinematic viscosity of the liquid. To understand this, we study the approach
of a liquid drop towards a solid surface in the presence of an intervening gas
layer. We take a numerical approach to solve the Navier–Stokes equations for
the liquid, coupled to the compressible lubrication equations for the gas, in
two dimensions. With this approach, we recover the experimentally captured
early time effect of liquid viscosity on the drop impact, but our results show
that lift-off time and liquid kinematic viscosity have a more complex
dependence than the square root scaling relationship. We also predict the
effect of interfacial tension at the liquid–gas interface on the drop impact,
showing that it mediates the lift-off behavior.

###### keywords:

drops, capillary flows, Navier–Stokes equations, computational methods

##  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
(1877), 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 several review articles
(Rein, 1993; Yarin, 2006; Josserand & Thoroddsen, 2016).

Physical parameters that are known to affect the dynamics of a drop upon
impact with a substrate include the initial size and velocity of the drop
(Pasandideh-Fard et al., 1996), viscosity of the liquid (Pasandideh-Fard et
al., 1996), pressure and viscosity of the surrounding medium (Xu et al.,
2005), interfacial tension between the liquid and its surrounding medium
(Rioboo et al., 2003; Pasandideh-Fard et al., 1996), and the physical
properties of the substrate (Range & Feuillebois, 1998; Rein, 1993; Yarin,
2006; Josserand & Thoroddsen, 2016). Experiments by Xu et al. (2005) 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. 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.
The series of theoretical investigations by Mandre, Mani and Brenner (Mandre
et al., 2009; Mani et al., 2010; Mandre & Brenner, 2012) considers an inviscid
liquid drop falling towards a solid surface in the presence of an ideal gas,
and predicts 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 modeled using values of relevant field
variables exclusively at the liquid–gas interface, 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).

Experimental studies confirm 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, and show measurements of the air layer
underneath the impacting drop (Kolinski et al., 2012; van der Veen et al.,
2012; de Ruiter et al., 2012). However, subsequent experimental work by
Kolinski et al. (2014b) 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., 2009; Mani et al., 2010;
Mandre & Brenner, 2012).

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 drop. Previous computational
work on drop impact has been concerned with the dynamics of spreading (Eggers
et al., 2010) and splashing (Josserand & Zaleski, 2003; Boelens & de Pablo,
2018), whereas we focus on the dynamics of the drop in its approach towards a
solid surface, prior to contact.

We solve the incompressible Navier–Stokes equations in the drop, using a
modern implementation of Chorin’s projection method (Chorin, 1967, 1968) that
incorporates improvements from Almgren, Bell, Collela and coworkers (Bell et
al., 1989; Colella, 1990; Almgren et al., 1996). 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 experiment. We validate our model using
theoretical results (Mandre et al., 2009; Mani et al., 2010; Mandre & Brenner,
2012). 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.
(2014b). 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., 2009) as well as
experimental results (Kolinski et al., 2014b). 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 previous theoretical work (Mandre et
al., 2009; Mani et al., 2010; Mandre & Brenner, 2012). 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.

Figure 1: 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
viscosity is important. (c) Control volume and grid for discretization
constructed in the fluid domain.

###  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
modeling 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.

Symbol | Quantity | Order of magnitude | Units
---|---|---|---
$P_{0}$ | Initial gas pressure | $\mathcal{O}(10^{5})$ | $\mathrm{Pa}$
$R$ | Initial drop radius | $\mathcal{O}(10^{-3})$ | $\mathrm{m}$
$V_{0}$ | Initial drop velocity | $\mathcal{O}(10^{-1}\text{--}10^{0})$ | $\mathrm{m/s}$
$\gamma$ | Ratio of specific heat capacities, gas | $\mathcal{O}(1)$ | –
$\rho_{g}$ | Density of gas | $\mathcal{O}(10^{0})$ | $\mathrm{kg/m^{3}}$
$\rho_{l}$ | Density of liquid | $\mathcal{O}(10^{3})$ | $\mathrm{kg/m^{3}}$
$\nu_{g}$ | Kinematic viscosity of gas | $\mathcal{O}(10^{-5})$ | $\mathrm{m^{2}/s}$
$\nu_{l}$ | Kinematic viscosity of liquid | $\mathcal{O}(10^{-6}\text{--}10^{-4})$ | $\mathrm{m^{2}/s}$
Table 1: 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.

A theoretical consideration of our physical set-up for the parameter regime
specified in Table 1. In this parameter regime, the potential flow theory
argues that the relevant Reynolds number, based on the initial velocity
$V_{0}$ and drop radius $R$, is $\mathcal{O}\left(10^{2}\right)$, 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 build up 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., 2012; van der
Veen et al., 2012; de Ruiter et al., 2012) 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. (2014b) 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., 2009; Mani et al., 2010; Mandre & Brenner, 2012). 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$ 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. (2014b) observe that $\tau_{c}\propto{\nu_{l}}^{1/2}$.

As a consequence of the discrepancy between the observations of Kolinski et
al. (2014b), 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. (2014b), 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 co-ordinates, and those of the gas layer
in terms of three-dimensional cylindrical polar co-ordinates. In this co-
ordinate system, from the potential flow theory (Mandre et al., 2009; Mani et
al., 2010; Mandre & Brenner, 2012) and the experimental observations of air
layer profiles under an impacting drop (Kolinski et al., 2012; de Ruiter et
al., 2012; van der Veen et al., 2012; Kolinski et al., 2014b, a, 2019), we use
the simplifying approximation of no azimuthal variation in the flow field of
the drop or the gas.

As in the potential flow model, 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 co-ordinates $(x,y)$. The initial height
$h(x,t=T_{0})$ of the liquid–gas interface is described by a parabolic profile
as

$h(x,t=T_{0})=H_{0}+\frac{x^{2}}{2R},$ (1)

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

$\nabla\cdot\mathbf{u}_{l}=0,$ (2)

and the momentum equation for the liquid is

$\mathbf{u}_{l,t}+\mathbf{u}_{l}\cdot\mathbf{\nabla}\mathbf{u}_{l}=-\frac{\mathbf{\nabla}p_{l}}{\rho_{l}}+\nu_{l}\nabla^{2}\mathbf{u}_{l},$
(3)

where $\mathbf{u}_{l}=\left(u_{x},u_{y}\right)$, $p_{l}$ and $\rho_{l}$ are
the velocity field and pressure field in the liquid, respectively.

Based on the potential flow theory (Mandre et al., 2009; Mani et al., 2010;
Mandre & Brenner, 2012) and subsequent experimental observations (de Ruiter et
al., 2012; van der Veen et al., 2012; Kolinski et al., 2014b), 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,

$12\mu_{g}\left(\rho_{g}h\right)_{t}=\left(\rho_{g}h^{3}p_{g,x}\right)_{x},$
(4)

and the equation of state is described by an adiabatic assumption,

$\frac{p_{g}}{P_{0}}=\left(\frac{\rho_{g}}{\rho_{0}}\right)^{\gamma}.$ (5)

The coupling of flows in the liquid and gas at their interface is done using
Gibbs’ condition for pressure,

$p_{l}=p_{g}+\sigma h_{xx},$ (6)

and equality of shear stress in the liquid and gas at the liquid–gas interface
is

$\tau_{l}(x,h)=\tau_{g}(x,h).$ (7)

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
$\left(u,v_{x}\right)=(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 Section 2.1 and the mathematical
model, with initial and boundary conditions, in Section 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 sections 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 at the interface, close the the middle, 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 (Section 2.4). With the use of symmetry in the domain, as described in
Section 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}\mathbf{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}\mathbf{u}_{l}$ makes a
negligible contribution to equation (3). A schematic of the region where
viscosity is important in the deformed drop is shown by yellow in Figure 1(b)
and (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
equations (6) and (7). Equation (7) is enforced through a condition on the
gradient of the horizontal velocity,

$u_{l,y}(x,0)=\frac{\mu_{g}}{\mu_{l}}u_{g,y}(x,0).$ (8)

####  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, equation (3), by noting that the viscous term,
$\nu_{l}\nabla^{2}\mathbf{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, ${\bf u}_{l}\cdot\nabla{\bf u}_{l}$ in equation (3), will
also be small. With these approximations, we obtain a simplified momentum
equation for the liquid at these boundaries,

$\mathbf{u}_{l,t}=-\frac{\mathbf{\nabla}p_{l}}{\rho_{l}}.$ (9)

We use this simplified momentum equation to obtain velocity boundary
conditions at the interior boundaries. Taking the divergence of equation (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,

$p_{l}(x,y)=\frac{1}{\pi}\int_{-\infty}^{\infty}p_{l}(s,0)\frac{y}{(x-s)^{2}+y^{2}}\
\mathrm{d}s.$ (10)

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 equation
(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\leq t\leq t^{\text{end}}$, for a simulation domain of $0\leq x\leq L$ and
$0\leq y\leq\beta L$, where $\beta$ is the aspect ratio, for the liquid, and a
domain of $0\leq x\leq 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., 2009). 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

$\epsilon=\frac{P_{0}R\textit{St}^{4/3}}{\mu_{g}V},$ (11)

which is the ratio of the initial gas pressure to the pressure that would have
built up in the gas were treated as incompressible.

The results of Kolinski et al. (2014b), 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

$H^{*}=R\textit{St}^{2/3}$ (12)

where $\textit{St}=\mu_{g}/(\rho_{l}VR)$ is defined as the inverse of a
modified Stokes number.111The inverse of the modified Stokes number,
$\textit{St}=\mu_{g}/(\rho_{l}VR)$, has properties of two different media,
where the usual Stokes number is defined for a single medium. This scaling
result is strictly valid for inviscid flow, $\nu_{l}=0$,
$\Rey\rightarrow{}\infty$, and no surface tension at the liquid–air interface,
$\sigma=0$, $\textit{We}\rightarrow{}\infty$. While our simulations
incorporate liquid viscosity and surface tension, equation (12) still provides
a good estimate for the height at which the stagnation point forms, and is
useful measure of the physical scales of interest. We therefore set the
initial height to be a multiple of $H^{*}$, so that

$H_{0}=\tilde{H}_{0}H^{*}=\tilde{H}_{0}R\textit{St}^{2/3}$ (13)

for a dimensionless constant $\tilde{H}_{0}$. Similarly, we set
$t_{\text{end}}=\tilde{t}_{\text{end}}R\textit{St}^{2/3}/V$, and based on the
parabolic shape of the initial profile we choose
$L=\tilde{L}R\textit{St}^{1/3}$, where $\tilde{L}$ and
$\tilde{t}_{\text{end}}$ are dimensionless constants.

##  3 Numerical implementation

In this section, we describe the numerical implementation of the mathematical
model described in Section 2.2. Throughout the liquid domain, we compute and
record the liquid velocity $\mathbf{u}_{l}$ and pressure $p_{l}$. In the gas
layer, compute and track the height $h$ and pressure $p_{g}$. A reader may
skip the remainder of section without loss of continuity.

The flow-field variables $\phi$ are computed at discretized time steps, with
simulation time $t=n\Delta t$, where $\Delta t$ is the time step
discretization used in the simulation, and $n$ is an integer counter. Given
values of field variables $\phi^{(n)}=\phi(t=n\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, $\mathbf{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,
1967, 1968), with further developments in the computational techniques behind
the fluid simulation method (Bell et al., 1989; Almgren et al., 1998). In
particular, the work of Almgren et al. (1996) introduces the approximate
projection method discretized using the finite-element method (FEM). The
numerical methods are described in detail by Sussman et al. (1999), Yu et al.
(2003, 2007), and Rycroft et al. (2020). Here, we sketch the main features of
these methods.

We solve for $\mathbf{u}_{l}$ and pressure $p_{l}$ in the liquid domain by
solving the momentum equation for the liquid, equation (3), subject to
incompressibility of the liquid, as specified in equation (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 $\Delta x$ by $\Delta y$, as shown in Figure 2. The
velocities $\mathbf{u}^{n}$ are located in the cell-centers, 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 $\mathbf{u}^{*}$ according to

$\frac{{\bf u}^{*}-{\bf u}^{n}}{\Delta t}=-\left[{\bf u}\cdot\nabla{\bf
u}\right]^{n+1/2}+\frac{\nu}{2}\left(\nabla^{2}{\bf u}^{n}+\nabla^{2}{\bf
u}^{*}\right).$ (14)

Here, the viscous stress terms $\nu\nabla^{2}\mathbf{u}$ and
$\nu\nabla^{2}\mathbf{u}^{*}$ evaluated using a standard five-point finite-
difference stencil. Due to the presence of the ${\bf u}^{*}$ on the right hand
side of equation (14) it must be solved implicitly. We use a multigrid method
to solve the resulting linear system.

The advective term $[\mathbf{u}\cdot\nabla\mathbf{u}]^{n+1/2}$ is evaluated at
the half-timestep $n+1/2$ using the second-order Godunov upwinding scheme of
Colella (1990). 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 center, which is done using the fourth-order monotonicity-limited scheme
of Colella (1985). 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 et al.,
1991; Sussman et al., 1999). After this, the advective term can be accurately
evaluated using centered finite differences of the upwinded edge-based
velocities (Sussman et al., 1999; Yu et al., 2003, 2007; Rycroft et al.,
2020).

The velocity at step $n+1$ can be calculated from the intermediate velocity by
evaluating the pressure gradient term,

$\frac{{\bf u}^{n+1}-{\bf u}^{*}}{\Delta t}=-\frac{1}{\rho}{\nabla p^{n+1}},$
(15)

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 equation (15) and enforce that
$\nabla\cdot\mathbf{u}^{n+1}=0$, according to equation (2). Hence,

$\nabla\cdot\left[\Delta t\left\\{\frac{1}{\rho}\nabla
p^{n+1}\right\\}\right]=\nabla\cdot{\bf u}^{*},$ (16)

which is a Poisson equation for the pressure that can be solved. Once
$p^{n+1}$ is known, equation (15) can be used to evaluate $\mathbf{u}^{n+1}$
and complete the timestep from $n$ to $n+1$. We solve equation (16) using the
FEM discretization of Almgren et al. (1996). Here, each pressure value at a
cell corner has a corresponding bilinear hat function over the four
neighboring grid cells (Yu et al., 2003, 2007; Rycroft et al., 2020). Once the
pressure is computed, the gradient in equation (15) is evaluated using
centered finite differences of the pressure field values.

Figure 2: Discretization of the field variables. A: The two-dimensional
computational grid for the liquid. B: The one-dimensional computational grid
for the gas.

###  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 equation (14).

The velocity boundary conditions at the top and right boundaries of the liquid
domain, corresponding to the interior of the drop, and described in Section
2.3.3, are obtained through substituting the expression for pressure given by
equation (10) in the expression for velocity at these boundaries, equation
(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 $\mathbf{u}$ in the
liquid.

As shown in Appendix A.1, substituting the equation of state (5) into the
lubrication equation (4) yields

$12\frac{\mu}{\gamma}hp_{t}+12\mu
h_{t}p=\frac{1}{\gamma}h^{3}p_{x}p_{x}+3h^{2}h_{x}pp_{x}+h^{3}pp_{xx},$ (17)

where the $g$ subscript on the pressure has been dropped. We solve for this
equation as follows:

  1. 1.

$p$ is known at the time step $t^{(n)}=n\ \Delta t$

  2. 2.

$h$, $h_{t}$, and therefore $h_{x}$ are known at time step $t^{(n)}=n\ \Delta
t$

  3. 3.

solve for $p$ at time step $t^{(n+1)}=(n+1)\ \Delta t$

  4. 4.

use $p$ as a boundary condition for the flow in the liquid, to obtain
$v=h_{t}$ at time step $t^{(n+1)}$

To numerically solve equation (17) we discretize the spatial derivatives using
second-order centered 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

$12\frac{\mu}{\gamma}\bar{h}\left(\frac{p^{n+1}_{i}-p^{n}_{i}}{\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\Delta
x}\right)^{2}\\\
+3\bar{h}^{2}\bar{h}_{x}p^{n+1}_{i}\left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2\Delta
x}\right)+\bar{h}^{3}p^{n+1}_{i}\left(\frac{p^{n+1}_{i+1}-2p^{n+1}_{i}+p^{n+1}_{i-1}}{\Delta
x^{2}}\right),$ (18)

that can be solved for the gas pressure $p^{n+1}_{i}$. In equation (18), 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 centered
differencing, so that

$\bar{h}=\frac{h_{i}+h_{i+1}}{2},\qquad\bar{h}_{t}=\frac{h_{t,i}+h_{t,i+1}}{2},\qquad\bar{h}_{x}=\frac{h_{i+1}-h_{i}}{\Delta
x}.$ (19)

Due to the products of pressure terms on the right hand side of equation (18),
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 $\Delta
x=\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 (6) is small, there is
no timestep restriction in the simulation that scales like $\Delta x^{2}$
(Heath, 2002). We therefore choose a candidate timestep to satisfy $\Delta
t=\zeta\Delta x$ for a dimensionless constant $\zeta$, based on satisfying
Courant–Friedrichs–Lewy conditions (Courant et al., 1967) 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_{\text{end}}/N_{f}$. In general, an
integer multiple of candidate timesteps will not exactly match $t_{f}$, so
that $c\Delta t<t_{f}<(c+1)\Delta t$. Because of this, the actual timestep is
slightly adjusted to $\Delta t^{\prime}=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.

| Quantity | Value | Units
---|---|---|---
$P_{0}$ | Initial gas pressure | $10^{5}$ | $\mathrm{Pa}$
$R$ | Initial drop radius | $1.5\times 10^{-3}$ | $\mathrm{m}$
$V$ | Initial drop velocity | $0.45$ | $\mathrm{m/s}$
$\gamma$ | Ratio of specific heat capacities, gas | $1.4$ | –
$\sigma$ | Interfacial surface tension | $0.072$ | $\mathrm{N/m}$
$\rho_{g}$ | Density, gas |  $1.2754$ (density of air at 20 ∘C) | $\mathrm{kg/m^{3}}$
$\rho_{l}$ | Density, liquid | $997.96$ | $\mathrm{kg/m^{3}}$
$\nu_{g}$ | Kinematic viscosity, gas | $1.506\times 10^{-5}$ | $\mathrm{m^{2}/s}$
Table 2: 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.

##  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 Section 2.4 for the height of the stagnation point $H^{*}$ (Mandre et al.,
2009; Mani et al., 2010; Mandre & Brenner, 2012). We use the baseline physical
parameters given in Table 2. Since we only want to resolve initial
deceleration of the drop, and not the formation of the thin gas layer, we use
the computational parameters in column 3 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^{*})$.

| Quantity |  Initial dynamics, all values of $\nu_{l}$ |  Lift-off, $\nu_{l}\leq 20\,\text{mm}^{2}\\!/\text{s}$ |  Lift-off, $\nu_{l}>20\,\text{mm}^{2}\\!/\text{s}$
---|---|---|---|---
$(M,N)$ | Grid dimensions | $(2048,256)$ | $(5120,960)$ | $(5120,960)$
$\beta$ | Aspect ratio | $\nicefrac{{1}}{{8}}$ | $\nicefrac{{16}}{{3}}$ | $\nicefrac{{16}}{{3}}$
$\zeta$ | Timestep multiplier | $8\times 10^{-3}$ | $8\times 10^{-3}$ | $8\times 10^{-3}$
$\tilde{L}$ | Rescaled domain width | $30$ | $30$ | $45$
$\tilde{H}_{0}$ | Rescaled initial drop height | $15$ | $15$ | $15$
$\tilde{t}_{\text{end}}$ | Rescaled duration | $50$ | $50$ | $100$
$N_{f}$ | Number of frames | $250$ | $250$ | $250$
Table 3: 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 Section 4.1. The second two sets are for the
lift-off calculations in all other sections. The procedure for connecting the
rescaled variables $\tilde{L}$, $\tilde{H}_{0}$, and $\tilde{t}_{\text{end}}$
to physical values is described in Section 2.4.

As described in Section 2.4, at low initial velocities the flow in the gas
layer can be treated as incompressible, and the scaling result
$H^{*}=R\text{St}^{2/3}$ can be derived. Mani et al. (2010) extend this
analysis to look at higher initial velocities, where compressibility of the
gas becomes important. This happens at a height of

$H_{*}=R(\mu_{g}V/RP_{0})^{1/2},$ (20)

where the subscript ‘$*$’ is used to distinguish from the stagnation height.
By modeling the subsequent drop deceleration below $H_{*}$ assuming gas
compressibility, Mani et al. (2010) derive the result

$H^{*}_{\text{comp}}=R\textit{St}^{2/3}\epsilon^{(2-\gamma)/(2\gamma-1)}$ (21)

for the stagnation height.

We first performed sequences of simulations using low initial velocites over
the range from $V=0.15\,\text{m}/\text{s}$ to $V=1.35\,\text{m}/\text{s}$,
using four different liquid viscosities from
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ to
$\nu_{l}=300\,\text{mm}^{2}\\!/\text{s}$. 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\textit{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\textit{St}^{2/3})$ is approximately constant and is in agreement
with equation (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. At small
$\epsilon^{-1}$ the values of $H^{*}/(R\textit{St}^{2/3})$ are slightly lower
for the cases of larger $\nu_{l}$ and $\sigma$. This hints at the limits of
the potential flow theory predictions (Mandre et al., 2009; Mani et al.,
2010).

We also performed two sequences of simulations with higher initial velocities
from $V=1\,\text{m}/\text{s}$ to $V=18\,\text{m}/\text{s}$ using
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$, to examine the compressible regime.
However, we note from equation (13) that our usual choice for initial height
scales according to $V^{-2/3}$ whereas from equation (20), 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}$ in equation (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 equation
(21).

Figure 3: (a) Plot of rescaled initial stagnation height
$H^{*}/(R\textit{St}^{2/3})$ as a function of the inverse of the dimensionless
compressibility parameter $\epsilon=P_{0}R\textit{St}^{4/3}/\mu_{g}V$, for a
range of different liquid viscosities $\nu_{l}$, surface tensions $\sigma$,
and gas parameters $\gamma$. When not otherwise stated, baseline parameters
from Table 2 are used. By default, the drop starts from a height $H_{0}$
scaled according to equation (13). For
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$, 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}$ into equation (13) are
shown. (b) Zoomed-in plot of the same data, showing the region bounded by the
dotted gray rectangle in panel (a).

###  4.2 Effect of liquid viscosity on drop profile

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 columns 4 & 5 of Table 3. Furthermore, since
the deviation of the interface from the solid surface happens later for high
viscosities (Kolinski et al., 2014b), we use a larger domain and longer
duration when $\nu_{l}>20\,\text{mm}^{2}\\!/\text{s}$, as indicated in column
5 of Table 3. We begin by using the baseline initial drop velocity of
$V=0.45\,\text{m}/\text{s}$ to match the experiments of Kolinski et al.
(2014b).

Figure 4 shows the height profiles for three different simulations using
liquid viscosities of $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$,
$\nu_{l}=32\,\text{mm}^{2}\\!/\text{s}$, and
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$. Panels (a), (b), and (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\,\text{\textmu{}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.

Figure 4: Profiles of the height $h(x,t)$ of the gas layer at intervals spaced
$6.004\,\text{\textmu{}s}$ apart, corresponding to an integer multiple of the
frame output interval $t_{f}$ described in Section 3.4, for liquid viscosities
of (a) $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$, (b)
$\nu_{l}=32\,\text{mm}^{2}\\!/\text{s}$, and (c)
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$. Panels (d), (e), and (f) show the
same data as (a), (b), and (c), respectively, but with a smaller range of $h$
to highlight the lift-off behavior. 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 with $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$. At
$t=24.02\,\text{\textmu{}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\,\text{\textmu{}m}$. By $t=72.05\,\text{\textmu{}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\,\text{\textmu{}s}$.

Figure 5: Snapshots of pressure $p$ (left) and vorticity
$\omega=[\nabla\times\mathbf{u}]_{3}$ for a portion of the liquid domain for a
simulation with liquid viscosity $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$. The
field values exhibit large variations in scale and also switch sign, so a
nonlinear mapping $f(\alpha)=\frac{1}{\log 10}\sinh^{-1}\frac{\alpha}{2}$ is
used to create the color 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\,\text{\textmu{}s}$ and
$t=96.07\,\text{\textmu{}s}$ the position of the front (given by the local
minimum of the profiles in Figure 4(a,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., 2014b). 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. (2014b).

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, make 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: Zoomed-in plots of the height of the gas layer at intervals spaced
$1.2009\,\text{\textmu{}s}$ apart for liquid viscosity
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$. 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 gray dashed lines show
the profiles with Gaussian smoothing applied, and the gray circles show the
global minima of the smoothed lines.

Figure 6(a) shows a further zoomed-in plot of the height profiles during lift-
off for the case of $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$. 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 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.

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}$ from
Table 2, we also consider reduced surface tension where the tip becomes
sharper. In this case, there can be numerical difficulties with in identifying
the minimum due to per-gridpoint variations 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\,\text{\textmu{}m}$. This results in a minimal alteration in the leading
tip position (Figure 6(b)), but improves robustness when analyzing 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.

Figure 7: 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}$, for a range of different
liquid viscosities. 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. Figure 8: 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}$, for a range
of different liquid viscosities. 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.

Figure 7 shows the trajectories of the leading tip for a range of liquid
viscosities from $\nu_{l}=2.5\,\text{mm}^{2}\\!/\text{s}$ to
$\nu_{l}=160\,\text{mm}^{2}\\!/\text{s}$. The moment of lift-off is indicated
on each trajectory, by identifying the global minimum of each curve. Two
regimes are visible. For $\nu_{l}\lessapprox 20\,\text{mm}^{2}\\!/\text{s}$,
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\,\text{mm}^{2}\\!/\text{s}$, 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\,\text{mm}^{2}\\!/\text{s}$, 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\,\text{mm}^{2}\\!/\text{s}$ and
$\nu_{l}=25\,\text{mm}^{2}\\!/\text{s}$. In Figure 8 we show trajectories for
the case when the initial drop velocity is doubled to
$V=0.9\,\text{m}/\text{s}$. The scale of $h$ is smaller in the plot, but the
relative amplitude of the capillary waves is increased. As a result they have
a noticeable effect of lift-off times at lower viscosities, with a large
difference in lift-off time between $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$
and $\nu_{l}=13\,\text{mm}^{2}\\!/\text{s}$.

We now compare to the experimental finding of Kolinski et al. (2014b) 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. (2014b) were not able to
directly view $h_{0}$, since the drop can only be observed once it enters an
evanescent field of height $h_{\text{ev}}=1\,\text{\textmu{}m}$. Instead, they
measured the first position and time $(h^{\prime},r^{\prime},t^{\prime})$ when
the drop enters the evanescent field, and assume that the drop is undeformed
at this position, so that $h^{\prime}=h_{0}-Vt^{\prime}+r^{\prime
2}/(2R^{2})$. Hence $h_{0}$ can be calculated, and therefore the time origin
is $t^{\text{alt}}_{0}=h_{0}/V-t^{\prime}+r^{\prime 2}/(2R^{2}V)$.

Figure 9: Liftoff time $\tau$ as a function of liquid viscosity $\nu_{l}$, for
a range of initial drop velocities $V$. The experimental data is taken from
Figure 4(a) of the paper by Kolinski et al. (2014b) (which used
$V=0.45\text{~{}m}/\text{s}$). Three data points are omitted from the plot for
low $V$ and low $\nu_{l}$ due to physical difficulties with running the
simulation—see Section 4.3.

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.
(2014b) are also plotted. Even though the experiments of Kolinski et al.
(2014b) are performed in a three-dimensional axisymmetric configuration, there
is good quantitative agreement with the two-dimensional simulation results.

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}$, the data
for low impact velocities looks better fit to $\nu_{l}^{\eta}$, where there
are two different values of $\eta$ with a change at $\nu_{l}\approx
20\,\text{mm}^{2}\\!/\text{s}$. This is also consistent with the two different
types of behavior 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}^{\text{alt}}$ as the time origin; this resulted in small shift
upwards of the data (since $t_{0}^{\text{alt}}<t_{0}$ in all cases), but did
not affect the overall patterns. As a third approach, we defined a time origin
$t_{0}^{\text{dat}}$ directly from the data, by finding the best fit to the
model
$\tau_{\text{dat}}=t-t_{0}^{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$
for the three free parameters $(t_{0}^{\text{dat}},\gamma,\alpha)$. Here
$\nu_{\text{sc}}=20\,\text{mm}^{2}\\!/\text{s}$ is chosen as an arbitrary
viscosity scale. Specifically, for the data points $(\nu_{l,k},t_{k})$, we
minimized the residual

$S(t_{0}^{\text{dat}},\gamma,\alpha)=\frac{1}{2}\sum_{k}\left(\log\gamma+\alpha\log\frac{\nu_{l,k}}{\nu_{\text{sc}}}-\log(t_{k}-t_{0}^{\text{dat}})\right)^{2}.$
(22)

We used the L-BFGS-B algorithm for bound-constrained nonlinear optimization
(Byrd et al., 1995), and enforced $t_{0}^{\text{dat}}<0.999t_{\text{min}}$ and
$\alpha>0$, where $t_{\text{min}}=\min_{k}\\{t_{k}\\}$. We started the
minimization using multiple initial guesses with
$t_{0}^{\text{dat}}\in[0.7t_{\text{min}},0.9t_{\text{min}}]$ and
$\alpha\in[0.4,1.2]$. For each pair $(t_{0}^{\text{dat}},\alpha)$ the value of
$\gamma$ is set so that the bracketed expression in equation (22) is zero for
the $\nu_{l}=20\,\text{mm}^{2}\\!/\text{s}$ 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}$ 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 is more consistent with a linear scaling relationship as opposed to
the square root relationship in Figure 9. Panels (b), (c), and (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=\tfrac{1}{2}$ to
match the exponent of Kolinski et al. (2014b) then the average residual
increases to $\bar{S}=0.3855$ and the data points exhibit an upward curve 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
Appendix D 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.

Figure 10: (a) Liftoff time $\tau_{\text{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 minimizing the three parameter
residual function in equation (22). (b–d) best fits of the parameters
$t_{0}^{\text{dat}}-t_{0}$, $\gamma$, and $\alpha$ for different initial drop
velocities. Figure 11: Lift-off time as a function of liquid viscosity, for a
range of (a) different drop radii, and (b) different surface tension values.
The data point for
$(\nu_{l},R)=(2.5\,\text{mm}^{2}\\!/\text{s},0.75\,\text{mm})$ is omitted due
to physical difficulties with running the simulation (Section 4.3). The data
points for $\sigma=0.0072\,\text{N}/\text{m}$ and $\nu_{l}\geq
100\,\text{mm}^{2}\\!/\text{s}$ are omitted because lift-off does not occur
over the simulation duration.

Figure 11(a) shows lift-off times (relative to $t_{0}$) for 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.

###  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 Section 4.2, we calculated the lift-off
times for $\sigma=0.0072\,\text{N}/\text{m}$,
$\sigma=0.036\,\text{N}/\text{m}$ and the $\sigma=0.144\,\text{N}/\text{m}$,
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}$ lift-off is eliminated for large values of
$\nu_{l}$, with the leading tip continuing to decrease over the course of the
simulation.

Figure 12: Profiles of the height of the gas layer at intervals spaced
$6.004\,\text{\textmu{}s}$ apart for liquid viscosities of (a)
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$, (b)
$\nu_{l}=32\,\text{mm}^{2}\\!/\text{s}$, and (c)
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$ using zero surface tension. Panels
(d), (e), and (f) show the same data as (a), (b), and (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.  Figure 13: Zoomed-in plots of the height of the
gas layer at intervals spaced $1.2009\,\text{\textmu{}s}$ apart for liquid
viscosities of $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ using zero surface
tension. 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 gray dashed lines show the profiles with Gaussian smoothing applied,
and the gray circles show the global minima of the smoothed lines.

Figure 11(b) strongly suggests that surface tension is important in creating
lift-off. Reducing from $\sigma=0.036\,\text{N}/\text{m}$ to
$\sigma=0.0072\,\text{N}/\text{m}$ 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\,\text{mm}^{2}\\!/\text{s}$,
$\nu_{l}=32\,\text{mm}^{2}\\!/\text{s}$, and
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$ 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\,\text{mm}^{2}\\!/\text{s}$, indicating that 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 Section 4.2 is sufficient to extract smooth leading tip
trajectories.

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 behavior. To test this, we
compared simulations with the baseline parameters to those on a finer grids
(Appendix C). Unlike the case for finite surface tension, the simulations on a
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 different
liquid viscosities $\nu_{l}$. For $13\,\text{mm}^{2}\\!/\text{s}<\nu_{l}\leq
40\,\text{mm}^{2}\\!/\text{s}$ we switched to a larger computational grid of
size $8192\times 1536$, and for $\nu_{l}\leq 13\,\text{mm}^{2}\\!/\text{s}$ we
switched to a very large grid of size $12288\times 2304$. For $\nu_{l}\leq
44\,\text{mm}^{2}\\!/\text{s}$ we were not able to simulate on a large enough
grid to adequately resolve the numerical diffusion, and hence results 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.

Figure 14: 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 and zero surface tension.
Simulations with $13\,\text{mm}^{2}\\!/\text{s}<\nu_{l}\leq
40\,\text{mm}^{2}\\!/\text{s}$ use a grid of size $8192\times 1536$, and
simulations with $\nu_{l}\leq 13\,\text{mm}^{2}\\!/\text{s}$ use a grid of
size $12288\times 2304$.

We also found a regime where the effect of surface tension can qualitatively
affect the lift-off behavior. 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}$ in the regime of low initial velocity,
$V=0.3\,\text{m}/\text{s}$ and low viscosity,
$\nu_{l}=6.5\,\text{mm}^{2}\\!/\text{s}$. 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\,\text{\textmu{}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 behavior further. The sharp features
visible in Figure 15(a) (e.g. at
$(x,h)\approx(700\,\text{\textmu{}m},100\,\text{\textmu{}m})$) cause the
simulation to terminate early, since the Newton–Raphson iterations can no
longer be solved to an acceptable level of accuracy.

Figure 15: (a) Profiles of the height of the gas layer at intervals spaced
$11.80\,\text{\textmu{}s}$ apart for liquid viscosity of
$\nu_{l}=6.5\,\text{mm}^{2}\\!/\text{s}$ and initial drop velocity of
$V=0.3\,\text{m}/\text{s}$, showing the development of capillary waves. (b)
Zoomed-in plot of the same data. The global minima of the curves are plotted
when they indicate the leading tip.

##  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., 2009; Mani et al., 2010; Mandre & Brenner, 2012). 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. (2014b). We have therefore
demonstrated that with the reduced model described in Section 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. (2014b). 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.

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. (2009); Mandre & Brenner (2012) and Mani et al.
(2010) 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 behavior 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 mesh refinement (Berger &
Oliger, 1984; Guittet et al., 2015), where the liquid–gas interface and the
region around the leading tip could use a finer mesh. This could 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., 2003, 2007), but
the routines that involve the gas layer would require modifications. For
example, the kernel in equation (10) would need to be modified. While the
current simulation is already in good agreement with the experimental results
of Kolinski et al. (2014b), 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., 2017). 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., 2017; Li et al., 2017).
They have also examined the case of nano-rough surfaces (Langley et al.,
2018). 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., 2019).

[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.

[Author ORCID]C. H. Rycroft, https://orcid.org/0000-0003-4677-6990; S. Mishra,
https://orcid.org/0000-0002-9487-5531; S. M. Rubinstein,
https://orcid.org/0000-0002-2897-2579

##  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, equation (5), into the
lubrication equation for the pressure in the gas layer, equation (4), we get

$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}.$
(23)

Dividing both sides by ${\rho_{0}}/{p_{0}^{1/\gamma}}$,

$12\mu_{g}\left(p_{g}^{1/\gamma}h\right)_{t}=\left(p_{g}^{1/\gamma}h^{3}p_{g,x}\right)_{x}.$
(24)

Using the product rule to expand the brackets,

$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}3h^{2}h_{x}p_{g,x}+p_{g}^{1/\gamma}h^{3}p_{g,xx}.$
(25)

Dividing both sides by $p_{g}^{1/\gamma-1}$,

$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}3h^{2}h_{x}p_{g,x}+p_{g}h^{3}p_{g,xx}.$
(26)

Dropping the subscript $g$ for gas yields

$12\mu\left(\frac{1}{\gamma}p_{t}h+ph_{t}\right)=\frac{1}{\gamma}p_{x}h^{3}p_{x}+p3h^{2}h_{x}p_{x}+ph^{3}p_{xx},$
(27)

which can be rearranged to give equation (17).

###  A.2 Numerical solution of the gas layer pressure

We now describe how to solve equation (18) 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 (18) as a nonlinear system $F(P)=0$, where the $i$th
component of $F$ is given by the $(\text{L.H.S.}-\text{R.H.S.})$ of equation
(18) evaluated at the $i$th gridpoint. An improved estimate $P^{k+1}$ for the
pressure is given in terms of $P^{k}$ by

$J_{F}(P^{k})\left(P^{k+1}-P^{k}\right)=-F(P^{k}),$ (28)

where $J_{F}(P^{k})$ is the Jacobian of $F$, which has components

$J^{k+1}_{F,ij}=\frac{\partial F_{i}}{\partial P^{k+1}_{j}}.$ (29)

Since the finite-difference stencils in equation (18) only involve adjacent
grid points, $J_{F}$ is a tridiagonal system. It can be written as

$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),$ (30)

where the terms are given by

$\displaystyle a_{i}=\frac{\partial F}{\partial p_{i-1}^{n+1}}$
$\displaystyle=\frac{1}{\gamma}\frac{\bar{h}^{3}}{\Delta
x}\left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2\Delta
x}\right)-\frac{3\bar{h}^{2}\bar{h}_{x}}{2\Delta
x}p^{n+1}_{i}+\frac{\bar{h}^{3}}{\Delta x^{2}}{p^{n+1}_{i}},$ (31)
$\displaystyle b_{i}=\frac{\partial F}{\partial p_{i}^{n+1}}$
$\displaystyle=-12\frac{\mu}{\gamma}\frac{\bar{h}}{\Delta
t}(1)-12\mu\bar{h}_{t}-\left[\frac{3\bar{h}^{2}\bar{h}_{x}}{2\Delta
x}\left({p^{n+1}_{i+1}-p^{n+1}_{i-1}}\right)\right.$
$\displaystyle\phantom{=}\left.+\frac{\bar{h}^{3}}{\Delta
x^{2}}\left(p_{i+1}^{n+1}-2p_{i}^{n+1}+p_{i-1}^{n+1}\right)+\frac{\bar{h}^{3}}{\Delta
x^{2}}p_{i}^{n+1}(-2)\right],$ (32) $\displaystyle c_{i}=\frac{\partial
F}{\partial p_{i+1}^{n+1}}$
$\displaystyle=-\frac{1}{\gamma}\frac{\bar{h}^{3}}{\Delta
x}\left(\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2\Delta
x}\right)-\frac{3\bar{h}^{2}\bar{h}_{x}}{2\Delta
x}p^{n+1}_{i}-\frac{\bar{h}^{3}}{\Delta x^{2}}{p^{n+1}_{i}}.$ (33)

Solving equation (28) can be done efficiently using LAPACK’s tridiagonal
solver dgtsv (Anderson et al., 1999). 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.

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 (10) requires applying
Simpson’s rule along the $M$ gridpoints on the bottom edge. This must be done
for the top edge of length $M$, and the right edge of length $N$, requiring
$\mathcal{O}(M(M+N))$ work in total. This is 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 $\mathcal{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 linear scaling with gridpoints, we would expect 81 ms to
increase to $(81\text{~{}ms})\tfrac{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
$\mathcal{O}(M)$) decreases. In particular, more V-cycles are required to
handle the implicit viscosity linear system when $\nu_{l}$ is larger.

|  Initial dynamics, $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ |  Lift-off, $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ |  Lift-off, $\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$
---|---|---|---
Grid size | $2048\times 256$ | $5120\times 960$ | $5120\times 960$
Total WC time | 0.641 h | 20.6 h | 29.8 h
Total timesteps | 28500 | 71250 | 95000
WC time $/$ timestep | 81 ms | 1041 ms | 1130 ms
WC fraction on gas layer | 0.25% | 0.047% | 0.045%
WC fraction on BCs | 3.39% | 1.61% | 1.49%
WC fraction on MAC solve | 11.21% | 13.05% | 12.09%
WC fraction on FEM solve | 12.77% | 14.31% | 14.36%
WC fraction on viscous solve | 20.11% | 26.67% | 31.12%
Mean MAC V-cycles | 5 | 4.4 | 4.5
Mean FEM V-cycles | 5.6 | 4.9 | 5.3
Mean viscosity V-cycles | 3.7 | 4.3 | 5.5
Table 4: 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 Table 2 & 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 marker-and-cell
(MAC) linear system, and solving finite-element method (FEM) linear system,
are reported. Figure 16: Plots of the height profiles at
$t=102.07\,\text{\textmu{}s}$ for simulations with (a) liquid viscosity
$\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ and (b) liquid viscosity
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$. Baseline parameters are used
although the original domain size for both values of $\nu_{l}$ uses
$\tilde{L}=30$ and $\beta=\nicefrac{{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 $\Delta x=\Delta y$ the same.

##  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 affect 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 $\Delta x=\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 16
shows the height profiles in the thin gas layer for the four simulations, for
(a) $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ and (b)
$\nu_{l}=100\,\text{mm}^{2}\\!/\text{s}$. The vertical extensions give
visually indistinguishable curves, indicating that vertical dimension is
sufficiently large to resolve the viscous effects even for the larger value of
$\nu_{l}$. The horizontal extensions create a 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.

Figure 17: Plots of the height profiles spaced $6.004\,\text{\textmu{}s}$
apart for three different resolutions using the baseline parameters and liquid
viscosity $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$, with (a) surface tension
$\sigma=0.072\,\text{N}/\text{m}$, and (b) zero surface tension.

We also examined the sensitivity of the results to the numerical grid size.
Figure 17(a) shows height profiles for a simulation using the baseline
parameters with liquid viscosity $\nu_{l}=10\,\text{mm}^{2}\\!/\text{s}$ and
surface tension $\sigma=0.072\,\text{N}/\text{m}$. Simulations with higher
resolutions of $8192\times 1536$ and $12288\times 2304$ are also plotted,
resulting in visually indistinguishable results. Figure 17(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 Section 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 behavior of the leading tip
difficult to resolve for small values of $\nu_{l}$, requiring larger grid
sizes as presented in Section 4.3.

##  Appendix D Additional model fits for the lift-off time

In Section 4.2 of the paper, we examine how the lift-off time $t$ varies as a
function of liquid viscosity $\nu_{l}$. For a sequence of measurements
$(\nu_{l,k},t_{k})$ we fit the model
$\tau_{\text{dat}}=t-t_{0}^{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$
for the three parameters $(t_{0}^{\text{dat}},\gamma,\alpha)$. Here
$\nu_{\text{sc}}=20\,\text{mm}^{2}\\!/\text{s}$ is an arbitrary viscosity
scale. We minimized the residual

$S(t_{0}^{\text{dat}},\gamma,\alpha)=\frac{1}{2}\sum_{k}\left(\log\gamma+\alpha\log\frac{\nu_{l,k}}{\nu_{\text{sc}}}-\log(t_{k}-t_{0}^{\text{dat}})\right)^{2},$
(34)

using the L-BFGS-B algorithm for bound-constrained nonlinear optimization Byrd
et al. (1995), and we enforced $t_{0}^{\text{dat}}<0.999t_{\text{min}}$ and
$\alpha>0$, where $t_{\text{min}}=\min_{k}\\{t_{k}\\}$. Figure 18 shows the
results of this optimization for the six different initial drop velocities
$V$, along with the values of the fitted parameters. The mean value of the
residual is $\bar{S}=0.254$.

In the experiments by Kolinski et al. (2014b), the parameter $\alpha$ is
proposed to be $\nicefrac{{1}}{{2}}$. We therefore ran an additional
optimization where $\nu_{l}$ is constrained to be $\nicefrac{{1}}{{2}}$. The
mean residual increases to $\bar{S}=0.385$ and the model does not capture the
curve in the data points (Figure 19). We ran a further optimization where
$\alpha=1$ (Figure 20), which results in $\bar{S}=0.305$. We found that the
value of $\alpha$ that minimizes the mean residual is $\alpha=1.1181$. For
this case $\bar{S}=0.281$ and the optimization results shown in Figure 21.

Figure 22 shows an optimization when the time origin is fixed to the original
proposed definition of $t_{0}=H/V$, where $H$ is the drop height above the
surface. In Figure 23 the exponent is additionally constrained to
$\alpha=\nicefrac{{1}}{{2}}$. For each of the six optimizations, the values of
$S$ are reported in Table 5.

Constraints | $S_{V=0.3\text{\,m}/\text{s}}$ | $S_{V=0.45\text{\,m}/\text{s}}$ | $S_{V=0.6\text{\,m}/\text{s}}$ | $S_{V=0.75\text{\,m}/\text{s}}$ | $S_{V=0.9\text{\,m}/\text{s}}$ | $S_{V=1.05\text{\,m}/\text{s}}$ | $\bar{S}$
---|---|---|---|---|---|---|---
– | 0.019 | 0.260 | 0.415 | 0.356 | 0.264 | 0.207 | 0.2536
$\alpha=\nicefrac{{1}}{{2}}$ | 0.226 | 0.483 | 0.523 | 0.461 | 0.332 | 0.287 | 0.3855
$\alpha=1$ | 0.093 | 0.466 | 0.417 | 0.358 | 0.266 | 0.229 | 0.3049
$\alpha=1.1181$ | 0.039 | 0.345 | 0.446 | 0.372 | 0.273 | 0.208 | 0.2805
$t_{0}^{\text{dat}}=t_{0}$ | 0.218 | 0.520 | 0.572 | 0.494 | 0.351 | 0.314 | 0.4116
$t_{0}^{\text{dat}}=t_{0},\alpha=\nicefrac{{1}}{{2}}$ | 0.227 | 0.522 | 0.602 | 0.578 | 0.488 | 0.553 | 0.4952
Table 5: Values of the residual from equation (34) for the six different
optimizations considered. The first column shows the additional constraints
that are applied on the three parameters $(t_{0}^{\text{dat}},\gamma,\alpha)$.
The next six columns show the residual for the six different initial drop
velocities $V$. The final column shows the mean residual taken over the six
values of $V$. Figure 18: (a) Model fits for minimizing the residual $S$ for
the three parameters $(t_{0}^{\text{dat}},\gamma,\alpha)$. The straight lines
show the fitted model
$\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$. (b–d) values of
the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$. Figure 19: (a)
Model fits for minimizing the residual $S$ for the two parameters
$(t_{0}^{\text{dat}},\gamma)$ where $\alpha=\nicefrac{{1}}{{2}}$. The straight
lines show the fitted model
$\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$. (b–d) values of
the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$, with the dotted
line indicating the constrained parameter. Figure 20: (a) Model fits for
minimizing the residual $S$ for the two parameters
$(t_{0}^{\text{dat}},\gamma)$ where $\alpha=1$. The straight lines show the
fitted model $\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$.
(b–d) values of the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$,
with the dotted line indicating the constrained parameter. Figure 21: (a)
Model fits for minimizing the residual $S$ for the two parameters
$(t_{0}^{\text{dat}},\gamma)$ where $\alpha=1.1181$. The straight lines show
the fitted model $\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$.
(b–d) values of the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$,
with the dotted line indicating the constrained parameter. Figure 22: (a)
Model fits for minimizing the residual $S$ for the two parameters
$(\gamma,\alpha)$ where $t_{0}^{\text{dat}}=t_{0}$. The straight lines show
the fitted model $\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$.
(b–d) values of the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$,
with the dotted line indicating the constrained parameter. Figure 23: (a)
Model fits for minimizing the residual $S$ for the paramater $\gamma$, where
$t_{0}^{\text{dat}}=t_{0}$ and $\alpha=\nicefrac{{1}}{{2}}$. The straight
lines show the fitted model
$\tau_{\text{dat}}=\gamma(\nu_{l}/\nu_{\text{sc}})^{\alpha}$. (b–d) values of
the parameters $t_{0}^{\text{dat}}$, $\gamma$, and $\alpha$, with the dotted
lines indicating the constrained parameters.

## References

  * Almgren et al. (1998) Almgren, Ann S., Bell, John B., Colella, Phillip, Howell, Louis H. & Welcome, Michael L. 1998 A conservative adaptive projection method for the variable density incompressible Navier–Stokes equations. Journal of Computational Physics 142 (1), 1–46.
  * Almgren et al. (1996) Almgren, Ann S., Bell, John B. & Szymczak, William G. 1996 A numerical method for the incompressible Navier–Stokes equations based on an approximate projection. SIAM Journal on Scientific Computing 17 (2), 358–369.
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. & Sorensen, D. 1999 LAPACK Users’ Guide, 3rd edn. Philadelphia, PA: Society for Industrial and Applied Mathematics.
  * Bell et al. (1989) Bell, John B., Colella, Phillip & Glaz, Harland M. 1989 A second-order projection method for the incompressible Navier–Stokes equations. Journal of Computational Physics 85 (2), 257–283.
  * Bell et al. (1991) Bell, John B., Howell, Louis H. & Colella, Phillip 1991 An efficient second-order projection method for viscous incompressible flow. In 10th Computational Fluid Dynamics Conference.
  * Berger & Oliger (1984) Berger, Marsha J. & Oliger, Joseph 1984 Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics 53 (3), 484–512.
  * Boelens & de Pablo (2018) Boelens, Arnout M. P. & de Pablo, Juan J. 2018 Simulations of splashing high and low viscosity droplets. Physics of Fluids 30 (7), 072106.
  * Byrd et al. (1995) Byrd, Richard H., Lu, Peihuang, Nocedal, Jorge & Zhu, Ciyou 1995 A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing 16 (5), 1190–1208.
  * Chorin (1967) Chorin, Alexandre J. 1967 A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics 2 (1), 12–26.
  * Chorin (1968) Chorin, Alexandre J. 1968 Numerical solution of the Navier–Stokes equations. Mathematics of Computation 22 (104), 745–762.
  * Colella (1985) Colella, Phillip 1985 A direct Eulerian MUSCL scheme for gas dynamics. SIAM Journal on Scientific and Statistical Computing 6 (1), 104–107.
  * Colella (1990) Colella, Phillip 1990 Multidimensional upwind methods for hyperbolic conservation laws. Journal of Computational Physics 87 (1), 171–200.
  * Courant et al. (1967) Courant, Richard, Friedrichs, Kurt & Lewy, Hans 1967 On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11 (2), 215–234.
  * Eggers et al. (2010) Eggers, Jens, Fontelos, Marco A., Josserand, Christophe & Zaleski, Stéphane 2010 Drop dynamics after impact on a solid wall: Theory and simulations. Physics of Fluids 22 (6), 062101\.
  * Guittet et al. (2015) Guittet, Arthur, Theillard, Maxime & Gibou, Frédéric 2015 A stable projection method for the incompressible navier–stokes equations on arbitrary geometries and adaptive quad/octrees. Journal of Computational Physics 292, 215–238.
  * Heath (2002) Heath, Michael T. 2002 Scientific Computing: An Introductory Survery. McGraw-Hill.
  * Josserand & Thoroddsen (2016) Josserand, Christophe & Thoroddsen, Sigurdur T. 2016 Drop impact on a solid surface. Annual Review of Fluid Mechanics 48, 365–391.
  * Josserand & Zaleski (2003) Josserand, Christophe & Zaleski, Stéphane 2003 Droplet splashing on a thin liquid film. Physics of Fluids 15 (6), 1650–1657.
  * Kolinski et al. (2019) Kolinski, John M, Kaviani, Ramin, Hade, Dylan & Rubinstein, Shmuel M 2019 Surfing the capillary wave: wetting dynamics beneath an impacting drop. Physical Review Fluids 4 (12), 123605.
  * Kolinski et al. (2014a) Kolinski, John M., Mahadevan, L. & Rubinstein, Shmuel M. 2014a Drops can bounce from perfectly hydrophilic surfaces. EPL (Europhysics Letters) 108 (2), 24001\.
  * Kolinski et al. (2014b) Kolinski, John M, Mahadevan, L. & Rubinstein, Shmuel M. 2014b Lift-off instability during the impact of a drop on a solid surface. Physical Review Letters 112 (13), 134501.
  * Kolinski et al. (2012) Kolinski, John M., Rubinstein, Shmuel M., Mandre, Shreyas, Brenner, Michael P., Weitz, David A. & Mahadevan, L. 2012 Skating on a film of air: drops impacting on a surface. Physical Review Letters 108 (7), 074503.
  * Langley et al. (2017) Langley, Kenneth, Li, Er Qiang & Thoroddsen, Sigurdur T. 2017 Impact of ultra-viscous drops: air-film gliding and extreme wetting. Journal of Fluid Mechanics 813, 647–666.
  * Langley et al. (2018) Langley, Kenneth R., Li, Er Qiang, Vakarelski, Ivan U. & Thoroddsen, Sigurdur T. 2018 The air entrapment under a drop impacting on a nano-rough surface. Soft Matter 14, 7586–7596.
  * Li et al. (2017) Li, Er Qiang, Langley, Kenneth R., Tian, Yuan Si, Hicks, Peter D. & Thoroddsen, Sigurdur T. 2017 Double contact during drop impact on a solid under reduced air pressure. Physical Review Letters 119, 214502.
  * Mandre & Brenner (2012) Mandre, Shreyas & Brenner, Michael P. 2012 The mechanism of a splash on a dry solid surface. Journal of Fluid Mechanics 690, 148–172.
  * Mandre et al. (2009) Mandre, Shreyas, Mani, Madhav & Brenner, Michael P. 2009 Precursors to splashing of liquid droplets on a solid surface. Physical Review Letters 102 (13), 134502.
  * Mani et al. (2010) Mani, Madhav, Mandre, Shreyas & Brenner, Michael P. 2010 Events before droplet splashing on a solid surface. Journal of Fluid Mechanics 647, 163–185.
  * Pasandideh-Fard et al. (1996) Pasandideh-Fard, M., Qiao, Y. M., Chandra, Sanjeev & Mostaghimi, Javad 1996 Capillary effects during droplet impact on a solid surface. Physics of Fluids 8 (3), 650–659.
  * Range & Feuillebois (1998) Range, Kai & Feuillebois, François 1998 Influence of surface roughness on liquid drop impact. Journal of Colloid and Interface Science 203 (1), 16–30.
  * Rein (1993) Rein, Martin 1993 Phenomena of liquid drop impact on solid and liquid surfaces. Fluid Dynamics Research 12 (2), 61.
  * Rioboo et al. (2003) Rioboo, R., Bauthier, C., Conti, J., Voue, M. & De Coninck, J. 2003 Experimental investigation of splash and crown formation during single drop impact on wetted surfaces. Experiments in Fluids 35 (6), 648–652.
  * de Ruiter et al. (2012) de Ruiter, Jolet, Oh, Jung Min, van den Ende, Dirk & Mugele, Frieder 2012 Dynamics of collapse of air films in drop impact. Physical Review Letters 108 (7), 074505.
  * Rycroft et al. (2020) Rycroft, Chris H., Wu, Chen-Hung, Yu, Yue & Kamrin, Ken 2020 Reference map technique for incompressible fluid–structure interaction. Journal of Fluid Mechanics 898, A9.
  * Sussman et al. (1999) Sussman, Mark, Almgren, Ann S., Bell, John B., Colella, Phillip, Howell, Louis H. & Welcome, Michael L. 1999 An adaptive level set approach for incompressible two-phase flows. Journal of Computational Physics 148 (1), 81–124.
  * van der Veen et al. (2012) van der Veen, Roeland CA, Tran, Tuan, Lohse, Detlef & Sun, Chao 2012 Direct measurements of air layer profiles under impacting droplets using high-speed color interferometry. Physical Review E 85 (2), 026315.
  * Worthington (1877) Worthington, Arthur Mason 1877 XXVIII. On the forms assumed by drops of liquids falling vertically on a horizontal plate. Proceedings of the Royal Society of London 25 (171–178), 261–272.
  * Xu et al. (2005) Xu, Lei, Zhang, Wendy W. & Nagel, Sidney R. 2005 Drop splashing on a dry smooth surface. Physical Review Letters 94 (18), 184505.
  * Yarin (2006) Yarin, Alexander L. 2006 Drop impact dynamics: splashing, spreading, receding, bouncing…. Annu. Rev. Fluid Mech. 38, 159–192.
  * Yu et al. (2003) Yu, Jiun-Der, Sakai, Shinri & Sethian, James A 2003 A coupled level set projection method applied to ink jet simulation. Interfaces and Free Boundaries 5 (4), 459–482.
  * Yu et al. (2007) Yu, Jiun-Der, Sakai, Shinri & Sethian, James A. 2007 Two-phase viscoelastic jetting. Journal of Computational Physics 220 (2), 568–585.
  * Zhang et al. (2019) Zhang, Weiqun, Almgren, Ann, Beckner, Vince, Bell, John, Blaschke, Johannes, Chan, Cy, Day, Marcus, Friesen, Brian, Gott, Kevin, Graves, Daniel, Katz, Max P., Myers, Andrew, Nguyen, Tan, Nonaka, Andrew, Rosso, Michele, Williams, Samuel & Zingale, Michael 2019 AMReX: a framework for block-structured adaptive mesh refinement. Journal of Open Source Software 4 (37), 1370.