# Near-wall depletion and layering affect contact line friction of
multicomponent liquids

Michele Pellegrino  Berk Hess  Swedish e-Science Research center, Science for
Life Laboratory, Department of Applied Physics KTH, Stockholm 100 44, Sweden

###### Abstract

The main causes of energy dissipation in micro- and nano-scale wetting are
viscosity and liquid-solid friction localized in the three-phase contact line
region. Theoretical models predict the contact line friction coefficient to
correlate with the shear viscosity of the wetting fluid. Experiments conducted
to investigate such correlation have not singled out a unique scaling law
between the two coefficients. We perform Molecular Dynamics simulations of
liquid water-glycerol droplets wetting silica-like surfaces, aimed to
demystify the effect of viscosity on contact line friction. The viscosity of
the fluid is tuned by changing the relative mass fraction of glycerol in the
mixture and it is estimated both via equilibrium and non-equilibrium Molecular
Dynamics simulations. Contact line friction is measured directly by inspecting
the velocity of the moving contact line and the microscopic contact angle. It
is found that the scaling between contact line friction and viscosity is sub-
linear, contrary to the prediction of Molecular Kinetic Theory. The
disagreement is explained by accounting for the depletion of glycerol in the
near-wall region. A correction is proposed, based on multicomponent Molecular
Kinetic Theory and the definition of a re-scaled interfacial friction
coefficient.

Molecular Dynamics, Wetting, Droplets, Moving contact lines, Nanofluidics,
Viscosity

##  I Introduction

The dynamics of wetting and dewetting involves the motion of three-phase
contact lines, which is thus ubiquitous both in nature and engineering.
Deepening the physical understanding of moving contact lines is therefore
essential in order to improve industrial processes such as coating [1],
3D-printing [2] and boiling [3], as well as to provide insight for the
development of bio-inspired surfaces [4, 5]. The research on this topic
involves expertise from several communities and has produced a plethora of
mathematical models [6, 7, 8, 9, 10]. Despite differences between modeling
approaches, it is clear that viscosity and liquid-solid friction represent the
main sources of energy dissipation, and thus dominate the dynamics of contact
lines at sufficiently small length scales [11].

The effect of liquid-solid friction localized in the three-phases region, i.e.
contact line friction, has been studied extensively and it has been determined
to produce a deviation of the dynamic contact angle from its Young-Dupré
equilibrium value [12, 13, 14]. The prevalent first-principle explanation of
contact line friction emerges from the Molecular Kinetic Theory (MKT)
formulated by Blake and Haynes [15]: when a liquid-vapor interface impinges a
flat solid surface, the energy balance in the displacement of a three-phase
contact line involves on one hand the work of adhesion of liquid molecules on
the solid surface and on the other hand the attraction forces between liquid
molecules in the near-wall region. The latter component is in turn postulated
to correlate with the viscosity of the wetting liquid. It entails that the
contact line friction coefficient for a fluid-solid combination scales
linearly with the shear viscosity coefficient [16].

Recently, numerous experimental studies have probed the applicability and the
implications of MKT [17, 18, 19, 20, 21]. For instance, Duvivier et al.
observed the predicted linear scaling between viscosity and contact line
friction in the spreading dynamics of water-glycerol droplets over glass [22].
The scaling was confirmed by a following work which analyzed 20 separate
dynamic contact angle studies [23]. The experiments on silica involved
estimating the contact line friction coefficient directly from optical
measurements of the dynamic contact angle [24]; this approach might not fully
capture the bending effect of friction on the liquid-vapor interfaces, which
extends below the optical resolution limit [25, 26]. Carlson et al. pursued a
different approach and inferred the value of contact line friction by
reproducing the spreading rate of water-glycerol droplets on silica by
simulating Cahn–Hilliard-Navier-Stokes equations with ad hoc wetting boundary
conditions [12, 27]. In disagreement with MKT, the scaling between contact
line friction and viscosity that matched simulation results was found to be
sub-linear. More recently, Li et al. estimated the liquid-solid friction
coefficient indirectly from spreading rates and found no linear correlation
with viscosity [28]. The disagreement between experimental studies indicates
the limitations of optical microscopy in disentangling the effect of viscous
bending at the visible scale (apparent contact angle) from the effect of
contact line friction at the molecular scale (microscopic contact angle).

Recent advances in molecular dynamics simulations have provided methods to
directly observe and analyze the nanoscopic dynamics of moving contact lines
via numerical experiments [14, 29, 30, 31, 32]. Nevertheless, no molecular
dynamics study has been focused on studying the correlation between viscosity
and contact line friction. In this paper we intend to investigate the
aforementioned mismatch between experimental evidence by employing molecular
simulations to carefully estimate viscosity and contact line friction of
water-glycerol droplets wetting hydrophilic surfaces. Rheological studies have
found that the shear viscosity coefficient of water-glycerol mixtures changes
drastically depending on the mass fraction of glycerol, while equilibrium
interfacial and wetting properties are affected to a lesser extent [33, 34].
Therefore simulating water-glycerol mixtures allows us to isolate and control
the effect of viscosity, in addition to imitating a fluid-surface combination
utilized in real-world experiments.

The material properties of simulated water-glycerol solutions are found to be
generally consistent with experimental results and theoretical prediction.
This holds particularly true for the shear viscosity coefficient. Moreover,
upon simulating two-dimensional droplet spreading we observe a clear 1:1
relation between dynamic contact angle and contact line speed, which is a
signature of contact line friction. Albeit the contact line friction
coefficient is higher the higher the viscosity, the scaling law is sub-linear.
Upon inspecting the liquid density structure in the near-wall region it is
found that glycerol depletes the solid surface, allowing water to form an
adsorbed layer. This simple observation explains the sub-linear scaling
qualitatively. A quantitative explanation is obtained by considering
multicomponent MKT [35]; by simply correcting for the local near-wall density
of each fluid component, the linear scaling is recovered.

The paper is organized as follows. In section II we briefly introduce the
linear contact line mobility model based on contact line friction. In section
III we illustrate the types of Molecular Dynamics simulations performed in
this work and their goals. In section IV we present the results of molecular
simulations, in particular the calculation on viscosity and contact line
friction. In section V we discuss the implication of simulation results,
focusing on the molecular effects at the liquid-solid interface and the
comparison with experimental studies. Finally, we draw our conclusions in
section VI.

##  II Contact line friction model

Figure 1: (a) Schematic representation of the moving contact line of a two-
dimensional spreading droplet. The reported reference system is the one
adopted in molecular simulations. The opaque droplet represents the initial
configuration and the transparent one the final equilibrium configuration. (b)
Molecular details of the liquid-solid combination, reporting the location of
the partial charge $q$ of silica quadrupoles.

Contact line motion is modeled in a two-dimensional geometry. We also assume
that one of the two fluids is dense and the other is its vapor phase. The
contact point where the liquid-vapor impinges on the solid surface moves with
velocity $u_{cl}$ in the direction parallel to the solid surface (figure 1a).
The contact angle $\theta$ is allowed to deviate from its equilibrium value
$\theta_{0}$. Molecular Kinetic Theory expresses the velocity of the contact
point in terms of frequency and length of discrete molecular displacement
events. For small differences between the equilibrium and the dynamic contact
angle, or equivalently for sufficiently small contact line velocities, the
following linear mobility relation is usually adopted:

$\mu_{f}u_{cl}=\Big{(}\frac{k_{B}T}{\kappa^{0}\lambda^{3}}\Big{)}u_{cl}=\sigma\big{(}\cos\theta_{0}-\cos\theta\big{)}\;,$
(1)

$\mu_{f}$ being referred to as the contact line friction coefficient and
$\sigma$ being the liquid-vapor interface tension; $\lambda$ is the length of
molecular jumps and $\kappa^{0}$ is the equilibrium molecular jump rate.

The equilibrium jump rate can be reasonably assumed to depend on a per-
molecule jump activation free energy $\Delta g$ according to Arrhenius
equation: $k^{0}\propto\exp\\{-\Delta g/(k_{B}T)\\}$. The jump activation free
energy can be split into one component that quantifies the adsorption
(desorption) work on (respectively from) the solid surface $\Delta g_{s}$, and
one encompassing the interactions between neighboring liquid molecules $\Delta
g_{l}$ [16]. The latter is in turns related to the viscosity of the wetting
fluid $\eta$ according to Eyring theory: $\eta\propto\exp\\{\Delta
g_{l}/(k_{B}T)\\}/v$, being $v$ the volume of a fluid molecule [36]. Hence,
the equilibrium jump rate is decomposed as:

$\kappa^{0}=\kappa_{s}^{0}\kappa_{l}^{0}\propto\frac{\kappa_{s}^{0}}{\eta
v}\;,$ (2)

Note that $\kappa_{s}^{0}$ and $\kappa_{l}^{0}$ should not be interpreted as
jump rates themselves, but simply as characteristic frequencies determined by
the work of adsorption on the substrate on one hand and on the other hand by
viscosity. Combining equations 1 and 2, the contact line friction coefficient
is found to scale linearly with the shear viscosity coefficient:

$\mu_{f}\propto\frac{k_{B}T}{\kappa_{s}^{0}}\frac{v}{\lambda^{3}}\eta\;.$ (3)

The simplest way to obtain the contact line friction coefficient from
molecular simulations of dynamic wetting is to fit equation 1 to the contact
line speed and the dynamic contact angle, which are both extracted from the
location and the shape of the liquid-vapor interface. In order to avoid over-
fitting, the equilibrium contact angle and the surface tension can be obtained
independently from simple equilibrium simulations. The shear viscosity
coefficient can also be obtained independently, although it necessitates
particular care, as will be illustrated in the following section.

##  III Methods

###  III.1 Molecular Dynamics simulations

In this section we briefly describe the molecular simulations performed to
measure the material parameters of aqueous glycerol droplets and to study how
they wet silica-like substrates. All simulations are run with GROMACS 2022
[37]. Additional information on the force field and simulation parameters are
provided in A.

The SPC/E model is used to parameterize water molecules, while the force field
parameters of glycerol molecules are taken from OPLS-AA according to the study
by Jahn et al. [38, 39]. The solid substrate consists of a monolayer of silica
quadrupoles arranged in a hexagonal lattice with spacing 0.45 nm (figure 1b).
Oxygen atoms in each quadrupole molecule are partially charged with charge
$-q$, to which corresponds a charge $2q$ on silicon atoms. The lattice
monolayer geometry does not reproduce the topology of real silica, which can
either be a three-dimensional crystal or amorphous, but emulates the
electrostatic interactions that are fundamental to correctly reproduce wetting
of polar liquids such as water or glycerol. The partial charge $q$ can be
tuned to change the wettability of the surface, and thus obtain different
contact angles.

Liquid mixtures are produced at six different mass fractions of glycerol,
indicated with $\alpha_{g}=\\{0.0,0.2,0.4,0.6,0.8,1.0\\}$. Cubic liquid boxes
are obtained by randomly inserting a fixed number of water and glycerol
molecules, according to the mass fraction. Boxes are then equilibrated at
constant hydrostatic pressure and temperature first (NPT), and then at
constant volume and temperature (NVT), with $T=300$ K and $P=1$ bar.
Production runs are performed at constant volume.

Scheme | Label |  Goal |  |  $L^{0}_{x}\times L^{0}_{y}\times L^{0}_{z}$
---|---|---|---|---
|  |  |  | [nm3]
| SLAB |  Surface tension |  | $10\times 5\times 10$
| BOX I |  Viscosity (equilibrium) |  | $5\times 5\times 5$
| BOX II |  Viscosity (non-equilibrium) |  | $5\times 5\times 15$
| DROP I |  Equilibrium contact angle |  | $60\times 4.7\times 30$
| DROP II |  Contact line friction |  | $120\times 4.7\times 45$
| MENISCUS |  Depletion layer |  | $160\times 4.7\times 30.4$
Table 1: Type of systems simulated using Molecular Dynamics. The figures shown
in the “scheme” column outline the systems’ geometry in the xz plane. Dashed
lines represent periodic boundary conditions, while solid lines represent
purely repulsive walls preventing molecules from traveling across the
boundaries perpendicular to the $z$ direction. The “label” and “goal” columns
report the name with which the systems will be referred to in the article and
the observables or transport coefficients that are to be obtained from
simulations. The right-most column reports the reference sizes of simulation
boxes, i.e. the sizes of initialized systems. $L^{0}_{i}$ may differ from the
final edge length since equilibration at constant pressure necessarily changes
the boxes’ volume.

In order to fully characterize the bulk, interfacial and wetting properties of
the water-glycerol-silica combination, several different configurations of
bulk liquid, liquid-vapor interfaces and liquid-vapor-solid contact lines need
to be simulated. Table 1 offers a brief overview. Configurations of type SLAB
are used to measure the liquid-vapor surface tension from the difference
between lateral and normal components of pressure:

$\sigma=\frac{L_{z}}{2}\Big{\\{}P_{zz}-\frac{P_{xx}+P_{yy}}{2}\Big{\\}}\;.$
(4)

Configurations of type DROP I and DROP II are used to obtain respectively the
equilibrium and the dynamic contact angle. Droplets are prepared with a
contact angle $\theta>\theta_{0}$ and simulated until fully spread. The
initial contact angle is obtained by setting $q=-0.35$ $|e^{-}|$, which is
then switched to $q=-0.72$ $|e^{-}|$ to simulate spontaneous wetting.
Subsequently, the equilibrium contact angle is measured by fitting a circular
cap to the liquid-vapor interface. Dynamic contact angle and contact line
speed are measured by tracking the interface over time for 5 independent
replicas of the same droplet spreading simulation with initial conditions
sampled from a long equilibrium trajectory. Interface extraction is based on
maps of liquid density binned on-the-fly while simulating and stored on
structured rectangular grids; further details are presented in B.

Configurations of type MENISCUS are used to extract information relative to
near-wall liquid density layering. A meniscus configuration is more
advantageous compared to a droplet configuration when extracting equilibrium
interfacial information, as it allows simulating two independent liquid-solid
interfaces simultaneously.

###  III.2 Shear viscosity

The dispersion of values for the viscosity of SPC/E water reported in the
literature is substantial, ranging from 0.64 cP [40] to 0.91 cP [41]. This
hints towards an inherent difficulty in finding a robust procedure to
calculate viscosity of liquids using Molecular Dynamics. Given the objective
of our study, the measurement of shear viscosity deserves particular care. We
utilize both equilibrium and non-equilibrium simulations and validate the
former against the latter. Additional details on the techniques can be found
in C.

The equilibrium approach consists in simulating configurations of type BOX I
at thermodynamic equilibrium for a sufficiently long time. The shear viscosity
coefficient is then obtained by computing the integral of the off-diagonal
components of the pressure tensor, that via the following Einstein relation:

$\eta=\lim_{t\rightarrow\infty}\frac{1}{6}\frac{V}{k_{B}T}\frac{d}{dt}\sum_{i\neq
j}\expectationvalue{\Big{(}\int_{0}^{t}P_{ij}(\tau)d\tau\Big{)}^{2}}\;,$ (5)

where angular brackets represent the ensemble average over several replicas.

Conversely, the non-equilibrium approach consists in applying an external
forcing to systems of type BOX II and measuring the response of the fluid. The
forcing can have the form of a periodic acceleration field [40]:

$a_{x}(z)=\xi\cos(kz)\;,\quad k=\frac{2\pi}{L_{z}}\;.$ (6)

The fluid response is obtained directly by inspecting the flow profile in the
direction parallel to the acceleration field, which can be shown to relax
exponentially in time to:

$u_{x}(z)=\frac{\xi\rho}{\eta k^{2}}\cos(kz)\;.$ (7)

In contrast to the equilibrium approach, the viscosity could in principle
depend on the amplitude of the external acceleration. Therefore, a few values
of $\xi$ need to be tried to assess the consistency of results.

##  IV Simulation results

###  IV.1 Material parameters and transport coefficients

Table 2 reports the results obtained for the values of glycerol mass fraction
under study. Liquid density is obtained from simulations of BOX I systems at
constant temperature and hydrostatic pressure, that is after the box volume is
allowed to relax to equilibrium. Density is found to increase with the mass
fraction of glycerol, in quantitative agreement with literature ([42], see the
supplementary information).

Surface tension does not depend substantially on $\alpha_{g}$, implying a
neutral effect of glycerol molecules on the liquid-vapor interface. A
calculation of the accessible surface area of glycerol and water molecules at
the liquid-vapor interface corroborates the absence of crowding or depletion
of any of the two species. These results are partially is in contrast with
experimental studies, where the surface tension of pure glycerol is reported
to be about 10% lower than the one of pure water [33, 34]. We believe the
reason for this discrepancy is the choice of water model, which underestimates
the surface tension between pure water and vapour. We observe a slight
improvement using TIP4P/2005, a more accurate and computationally expensive
water model. In this case we observe a roughly 4% drop between the surface
tension of pure water and the one of pure glycerol. A sensitivity analysis on
the estimate of contact line friction shows that using TIP4P/2005 instead of
SPC/E would have a marginal effect on the scaling between contact line
friction and viscosity and even the discrepancy with experimental results is
not large enough to affect the conclusions of the study. The results of the
sensitivity analysis are reported in the supplementary information.

For sufficiently low glycerol concentrations the equilibrium contact angle
does not change significantly. Values for our guess of specific surface
energy, i.e. partial charge on silica oxygens interacting with liquid
molecules, produce contact angles that are in the range of the ones from the
experimental work of Carlson et al. [27] and Yada et al. [34] for treated
silica surfaces (e.g. silanized).

$\alpha_{g}$ | $\rho$ | $\sigma_{SPC/E}$ | $\sigma_{TIP4P/2005}$ | $\theta_{0}$
---|---|---|---|---
| [kg/m3] | [$10^{-2}$Pa$\cdot$m] | [$10^{-2}$Pa$\cdot$m] | [∘]
0.0 | 997.9$\pm$0.1 | 5.59$\pm$0.03 | 6.14$\pm$0.18 | 48.94$\pm$1.60
0.2 | 1046$\pm$0.2 | 5.56$\pm$0.08 | 6.12$\pm$0.23 | 47.62$\pm$1.77
0.4 | 1096$\pm$0.1 | 5.56$\pm$0.04 | 6.23$\pm$0.28 | 47.12$\pm$1.78
0.6 | 1150$\pm$0.5 | 5.60$\pm$0.12 | 6.05$\pm$0.33 | 51.44$\pm$0.44
0.8 | 1192$\pm$0.8 | 5.49$\pm$0.24 | 5.80$\pm$0.37 | 53.40$\pm$1.77
1.0 | 1238$\pm$1.2 | 5.89$\pm$0.29† | / | 76.06$\pm$0.47
Table 2: Density, surface tension (for SPC/E and TIP4P/2005) and equilibrium
contact angle over silica monolayers of liquid mixtures for the tested value
of glycerol mass fraction. Density is simply computed from systems of type BOX
I at equilibrium. The surface tension of pure glycerol ($\dagger$) is not
reported for TIP4P/2005 since it does not depend on the water model. Figure 2:
(a): Shear viscosity (in centiPoise) as a function of the non-equilibrium
forcing amplitude $\xi$ and the mass fraction of glycerol $\alpha_{g}$. The
markers with error bars indicate the results of non-equilibrium simulations,
while dashed horizontal lines indicate the results of equilibrium simulations,
i.e. Einstein relation; shaded areas represent the uncertainty of equilibrium
results. (b): Comparison between Molecular Dynamics results, experimental
results by Segur and Oberstar [43] and the best fit of the empirical model by
Cheng [44]; $\beta$ is defined in equation 8.

The calculation of viscosity via equilibrium and non-equilibrium simulations
is performed for each value of $\alpha_{g}$. The results are presented in
figure 2a. Although it was not possible to obtain sufficiently converged
results for large mass fractions of glycerol and small periodic perturbations,
equilibrium and non-equilibrium approaches provide consistent results for all
other combinations of $\alpha_{g}$ and $\xi$. In the following, the viscosity
values obtained using Einstein’s formula will be taken as reference, as they
don’t depend on any arbitrary forcing parameter. Furthermore, the simulation
results are compared with the empirical model by Cheng [44]:

$\frac{\log(\eta/\eta_{g})}{\log(\eta_{w}/\eta_{g})}=(1-\alpha_{g})+\beta(\alpha_{g})\;,\quad\beta(\alpha_{g})=\frac{ab\alpha_{g}(1-\alpha_{g})}{a\alpha_{g}+b(1-\alpha_{g})}\;,$
(8)

$a$ and $b$ being fitting parameters, and with the experimental measurements
by Segur and Oberstar [43]. Results are reported in figure 2b and show
quantitative agreement. This comparison constitutes evidence that the rheology
of simulated aqueous glycerol solutions reflects physical reality.

###  IV.2 Contact line friction

Contact line friction is estimated directly by fitting equation 1 to the post-
processed results of molecular simulations, with fitting parameters
$\\{\mu_{f},\theta_{0}\\}$ for each glycerol fraction. The equilibrium contact
angle can be obtained from the intercept of the linear fit between $u_{cl}$
and $\cos\theta$, or can be imposed to the value obtained separately from
equilibrium runs (table 2). In order to incorporate both sources of
information, a ‘lasso-like’ regression is performed by imposing a weak
constraint on the cosine of the contact angle, which corresponds to the
minimization of the following functional:

$\mathcal{L}(\mu_{f},\theta_{0})=||\mu_{f}\vec{u}_{cl}-\sigma(\cos\theta_{0}-\cos\vec{\theta})||_{2}+\zeta|\cos\theta_{0}-\cos\widehat{\theta}_{0}|\;,$
(9)

being $\widehat{\theta}_{0}$ the value of the equilibrium contact angle
reported in table 2. The notation $\vec{\cdot}$ indicates the list of values
obtained from molecular simulations. The contact line friction parameter
resulting from the minimization of $\mathcal{L}$ is found to be insensitive to
the penalty term; the penalty coefficient $\zeta$ is set to 0.1. Further
information on uncertainty quantification can be found in the supplementary
information.

Figure 3: (a): Contact line speed against the opposite of the cosine of the dynamic contact angle. Markers represent measurements sampled from MD simulations results, whereas solid lines represent the best fit of equation 1, with cost function given by equation 9. (b): Log-log plot of the contact line friction coefficients against the shear viscosity coefficients, scaled on the reference value for water. The dotted line represents the best fit of MD results, while the solid and the dashed lines correspond respectively to the scaling laws found by Duvivier et al. [22] and Carlson et al. [27]. $\alpha_{g}$ | $\eta$ | $\mu_{f}$
---|---|---
| [cP] | [cP]
0.0 | 0.69$\pm$0.01 | 3.77$\pm$0.1
0.2 | 1.33$\pm$0.02 | 5.43$\pm$0.1
0.4 | 2.51$\pm$0.04 | 9.50$\pm$0.2
0.6 | 7.10$\pm$0.2 | 18.5$\pm$0.4
0.8 | 45.7$\pm$1.1 | 58.8$\pm$1.3
Table 3: Bulk shear viscosity and contact line friction as a function of the
mass fraction of glycerol. Note that contact line friction has the same
physical units of viscosity.

Figure 3a shows the measurements of the contact line speed against the dynamic
contact angle, as well as the best fit of equation 9. Albeit much thermal
noise remains even after averaging over a few replicas, a bijective relation
between contact line friction and dynamic contact angle can clearly be
observed for all values of $\alpha_{g}$: this is a signature of contact line
friction. Furthermore, the linear MKT model appears suitable to capture the
correlation between contact angle and contact line speed. Deviations from
linear behavior due to advancing/receding asymmetry and high contact line
speed have been previously reported [30, 45]; however, given that only
advancing contact are simulated and that the spreading process is slow,
especially for large glycerol concentrations, a linear model does suffice.

The values of the estimated contact line friction coefficients are reported in
table 3. Note that no result is reported for $\alpha_{g}=1.0$. On one hand, it
was not possible to obtain a sufficiently wide range of contact line
velocities from systems of type DROPLET II, given that the increase in
viscosity caused the spreading rate to slow dramatically. On the other hand,
interface tracking proved too noisy on smaller systems of type DROPLET I.

As shown in the log-log plot 3b, contact line friction does not scale linearly
with viscosity: to the same relative increase in shear viscosity corresponds a
smaller increase in contact line friction. This can be also easily observed
from table 3. Figure 3b reports the scaling laws found by Duvivier et al. and
by Carlson et al. [22, 27]; the scaling obtained from molecular simulations
lies in between $\sim\eta$ and $\sim\eta^{0.5}$.

###  IV.3 Glycerol depletion of interfacial layers

Figure 4: (a): Semi-logarithmic plot of the density profiles of water and
glycerol as a function of the vertical coordinate, for the case
$\alpha_{g}=0.2$. The density of the two components is re-scaled on the total
average density in the bulk. The inset shows a snapshot of the near-wall
position and orientation of water and glycerol molecules. (b) Cumulative
density profile (equation 10) re-scaled on the total bulk density, showing
depletion of glycerol molecules in the liquid-solid interfacial layer. The
inset shows the section of the meniscus used to compute density profiles.

The last set of results we report in this section regards near-wall liquid
density layering. For the sake of brevity, we will refer to the results of
$\alpha_{g}=0.2$, but the same conclusions hold for all glycerol
concentrations. Figure 4a reports the normalized density profiles of water and
glycerol along the vertical direction. Layering can be observed close to solid
walls. Regarding water, 3 layers can be clearly spotted labeled 1, 2 and 3 in
the figure: a layer of adsorbed molecules (1), a layer of molecules that are
hydrogen-bonded with silica (2) and a third layer that bonds to molecules of
layer 2 (3, indicated by a dotted line). Similar density oscillations have
been observed experimentally for water wetting flat mica surfaces [46]. As for
glycerol, there appears to be no adsorbed layer; a layer with density slightly
larger than the bulk can be observed at around 0.35nm from silicon atoms (4).

The spatial resolution of the density profile is limited by the cell size of
the density binning grid, 0.2 nm. In order to detect whether the liquid-solid
interface is crowded/depleted by a molecular species we compute the relative
cumulative density along the vertical direction:

$\overline{\rho}(z)=\frac{\int_{0}^{z}\rho(\zeta)d\zeta}{\int_{0}^{L_{z}}\rho(\zeta)d\zeta}\;.$
(10)

The gap between cumulative profiles entails crowding of water in the liquid
layers close to the interface, or conversely depletion of glycerol. In the
next section we further discuss the implications of glycerol depletion and how
it may be incorporated into MKT.

##  V Discussion

###  V.1 Two-component contact line friction model

In real-world experiments, the viscosity of a wetting fluid can be tuned by
changing its chemical composition, e.g. by the addition of a solute; it is
reasonable to suppose that the effects of such a modification would not be
limited to macroscopic rheology, but would also extend to microscopic
kinetics. Furthermore, the molecular structure of a liquid in the proximity of
interfaces may differ substantially compared to the one in the bulk: as shown
in section IV.3, the depletion of one species can be noticeable even in the
case of a mixture of simple molecules wetting a flat surface.

MKT has been originally formulated for the case of a single-component liquid:
on the light of our results and of the considerations above it may be
insufficient to explain wetting dynamics of solutions. An attempt to
generalize MKT to multicomponent liquids has been proposed by Liang et al.
[35]. The model formulation stems from the decomposition of the total
irreversible dissipation work at the contact line into the contribution from
each molecular species. The contact line speed is therefore expressed as
linear combination of the wetting speed of each component; in the limit of the
molecular relaxation being much faster than the contact line speed, the
linearized mobility relation reads:

$u_{cl}=\sigma\sum_{i=1}^{n}\frac{\lambda^{3}_{i}\kappa^{0}_{i}}{\alpha_{i}^{*}k_{B}T}\big{(}\cos\theta_{0}-\cos\theta\big{)}\;,$
(11)

being $n$ the number of molecular species. The superscript $\cdot^{*}$ denotes
values that are localized at the liquid-wall interface, as opposed to the
bulk. In the case of water-glycerol mixtures $n=2$. Instead of denoting
glycerol and water with $i=\\{1,2\\}$ we will use the notation $\\{g,w\\}$;
moreover $\alpha_{g}^{*}<\alpha_{g}$, owing to the results of section IV.3.
The dependence on viscosity is ‘hidden’ in the equilibrium jump rate
$\kappa_{i}^{0}\propto\kappa_{s,i}^{0}/(v_{i}\eta^{*})$, where we highlight
that the specific molecular volume will be species-dependent and that
viscosity in the interfacial region may differ from the one in the bulk.

###  V.2 Correction to linear scaling

$\alpha_{g}$ | $\expectationvalue{\tau^{hb}_{w}}$ | $\expectationvalue{\tau^{hb}_{g}}$ | $\lambda^{3}/v_{w}$ | $\lambda^{3}/v_{g}$ | $\alpha_{g}^{*}$ | $\eta^{*}$
---|---|---|---|---|---|---
| [ps] | [ps] |  |  |  | [cP]
0.2 | 2.38$\pm$0.40 | 2.03$\pm$0.37 | 5.23$\pm$0.03 | 1.07$\pm$0.02 | 0.186$\pm$0.022 | 1.147$\pm$0.080
0.4 | 2.40$\pm$0.44 | 2.08$\pm$0.39 | 5.17$\pm$0.03 | 1.08$\pm$0.01 | 0.356$\pm$0.030 | 2.135$\pm$0.265
0.6 | 2.50$\pm$0.40 | 2.17$\pm$0.35 | 5.07$\pm$0.05 | 1.08$\pm$0.01 | 0.502$\pm$0.032 | 4.334$\pm$0.754
0.8 | 2.24$\pm$0.42 | 2.51$\pm$0.48 | 4.91$\pm$0.08 | 1.09$\pm$0.01 | 0.720$\pm$0.037 | 20.97$\pm$6.702
Table 4: Average hydrogen bond lifetime of water and glycerol molecules with
silica, alongside the parameters used to evaluate equation 12: the ratio
between the cube of the lattice spacing and the per-molecule volume of water
and glycerol, the local mass fraction of glycerol and the estimated local
viscosity.

We now discuss how to simplify equation 11 and obtain a correction to the
relation between contact line friction and viscosity. We consider the
depletion of glycerol in the interfacial region and the difference in excluded
volume between the two molecular species to be the primary factors to correct
for. First, we assume that the layering and depletion effects discussed in
section IV.3 propagate up to the contact line. This assumption is
substantiated by visual inspection of density maps (figure 5a). There are
strong indications that the main contribution to contact line displacement is
provided by the jumping motion of molecules that are temporarily hydrogen
bonded to the substrate [14, 45]. Therefore, we assume that
$\lambda_{w}\simeq\lambda_{g}\simeq d_{hex}=0.45$ nm, since the spacing of the
hexagonal silica quadrupole lattice determines where the hydrogen bonding
sites are located. We also assume that the equilibrium jump rate on the
substrate does not change substantially with the addition of glycerol. This
assumption is substantiated by the following observations: on the one hand,
molecules interacting with the surface spend most their time bonded with
surface molecules rather than in the transition state between jumps; on the
other hand, the typical lifetime of hydrogen bonds between water and silica
and between glycerol and silica does not depend substantially on glycerol
concentration.

Owing to the assumptions discussed above, we can simplify equation 11 and
obtain a new relation between a re-scaled friction coefficient and viscosity:

$\mu_{f}^{*}=\Big{[}\frac{\lambda^{3}}{\alpha_{g}^{*}v_{g}}+\frac{\lambda^{3}}{(1-\alpha_{g}^{*})v_{w}}\Big{]}\mu_{f}\propto\frac{k_{B}T}{\kappa^{0}_{s}}\eta^{*}$
(12)

We now introduce a further assumption: $\eta^{*}=\eta(\alpha_{g}^{*})$, i.e.
the change of interfacial viscosity is only due to the depletion of one
molecular species. This assumption is consistent with the local average
density model (LADM), whereby the local shear viscosity is identified with the
viscosity of the homogeneous fluid at the local average density, or in this
case at the average relative concentration [47]. We remark that this may be a
very blunt simplification. It is known that the addition of glycerol disrupt
the local structure of the hydrogen bonds network of bulk-like water [48],
causing viscosity to deviate from ideal solution theory in its dependence on
glycerol concentration. However, it is not clear how this effect would
influence viscosity in nanoconfinement, i.e. in the occurrence of depletion
and adsorption. Layering in nanoconfinement is known to greatly affect the
local configurational shear viscosity close to solid walls in single-component
liquids, requiring a more refined model then LADM [49]. Hence, assuming that
the depletion, i.e. local change of mass fraction, is the only factor
affecting local viscosity should be considered a ‘zero-order’ approximation.

We can utilize the fit of equation 8 to re-compute viscosity with mass
fraction values at the liquid-solid interface. The molecular volumes $v_{w}$
and $v_{g}$ are measured indirectly from the free volume and the volume
occupied by each species in systems of type BOX I, while $\alpha_{g}^{*}$ is
simply computed by counting molecules within 0.8 nm from the solid wall in
systems of type MENISCUS (dotted line of figure 4b). Provided the interfacial
viscosity, the local mass fraction and the molecular volumes, expression 12
can be evaluated. The values of the coefficients used to evaluate equation 12
are reported in table 4. Appendix D briefly explains how the coefficients are
obtained from molecular simulations.

Figure 5: (a) Zoom on the lower-right contact line of a MENISCUS system
showing the density map of water (purple) and glycerol (green). (b) Log-log
plot of the corrected contact line friction (scaled on the friction of water)
against the corrected interfacial viscosity (scaled on the viscosity of
water), showing linear scaling between the two coefficients. Error bars
correspond to the uncertainty (one standard deviation) in $\alpha_{g}^{*}$,
see supplementary information.

Values for the depletion-corrected friction coefficient and the interfacial
viscosity are presented in figure 5b. Linear scaling between the two
coefficients is re-obtained. A leading-order correction for depletion and
excluded volume effects is therefore enough to reconcile the results of
Molecular Dynamics simulation with Molecular Kinetic Theory.

###  V.3 Comparison with experimental studies

As mentioned in the introduction, some dynamic contact angle experiments have
found a linear correlation between contact line friction and viscosity. In
this section we address the source of the mismatch between experiments and
molecular simulations. It is very possible that in some cases contact line
friction would be mainly determined by liquid-liquid interactions rather than
liquid-surface indications. These are cases of viscous coalescence analogy
[11], occurring for example when a particularly viscous and homogeneous liquid
droplet spreads on a hydrophobic or slippery surface. In such cases it is not
unreasonable to expect that experiments would find a 1:1 correlation between
viscosity and friction. Regardless, the method to estimate the contact line
friction coefficient can by itself ‘screen’ molecular-scale effects with
viscous interface bending. If the coefficient is estimated by applying
equations 1 to a optically measured contact angle, then what the experiment is
effectively measuring is the proportionality between spreading speed and the
difference between the equilibrium and an apparent contact angle:

$\mu_{f}^{app}u_{cl}=\gamma\big{(}\cos\theta_{0}-\cos\theta_{app}\big{)}\;.$
(13)

By using Petrov-Petrov combined molecular-hydrodynamic expression for the
apparent contact angle of an advancing contact line, it is possible to roughly
estimate the relation between $\mu_{f}^{app}$ and $\mu_{f}$ [8]:

$\theta_{app}^{3}\simeq\theta^{3}(u_{cl};\mu_{f})+\frac{\eta
u_{cl}}{\gamma}\log\big(z_{obs}/l^{*}\big{missing})\;,$ (14)

where $z_{obs}$ is the distance from the contact line where the contact angle
is measured and $l^{*}$ is a microscopic length scale, which can be
interpreted as the lower-scale cutoff for the interface curvature. We can
substitute equation 1 into 14 for the relation $\theta(u_{cl})$ and 14 into 13
for the apparent contact angle. After linearizing for small wetting speeds we
are left with the following expression:

$\mu_{f}^{app}\simeq\mu_{f}+\frac{3\log\big(z_{obs}/l^{*}\big{missing})}{\theta_{0}}\eta\;.$
(15)

Two considerations can be drawn from equation 15: (i) If the contact angle is
measured sufficiently far from the contact line, then it is only reasonable to
observe a linear scaling between apparent contact line friction and viscosity.
(ii) The apparent contact line friction always over-estimates molecular
friction; the effect is exacerbated the further the measurement of the contact
angle is from the contact line, the more viscous the spreading liquid is and
the more hydrophilic the surface is.

This last point helps explaining why contact line friction coefficients
measured in molecular simulations are usually an order of magnitude smaller
than the ones reported in experimental studies using optical microscopy to
measure the dynamic contact angle. However, even when the dynamic contact
angle is treated as a hidden parameter and contact line friction is inferred
from spreading rates, the reported friction coefficient might still be
substantially larger the typical values obtained from simulations [12, 27, 28,
34]. We believe this is due to the unavoidable presence of surface defects
such as topographic roughness and chemical heterogeneity, which have a great
impact on contact line mobility [50, 51]. Surface defects are rarely studied
with molecular simulations of wetting dynamics since the scale of defects is
sometimes much larger than what can be simulated and because reproducing
realistic nanoscale defects often requires more advanced molecular modeling.

###  V.4 Implications for the wetting properties of mixtures and solutions

Heterogeneous complex liquids are known to exhibit an interesting rheology.
For example, fluids containing colloidal particles or polymers (nanofluids)
have been studied because of their peculiar wetting properties. Self-assembly
close to solid surfaces and migration away/towards high-shear regions can
either reduce or increase spreading friction [52]. In other cases, the
addition of polymers can affect bulk rheology (e.g. by inducing shear thinning
and viscoelastic behaviors) without affecting contact line friction [34]. Our
results indicate that non-trivial correlations between bulk flow and wetting
can also occur in the case of simple liquids when they are formed by multiple
molecular species.

On high-friction surfaces (i.e. in cases where $\mu_{f}\gg\eta$, where wetting
is thus substrate-dominated [11]), the preferential affinity of substrate
molecules for one of the liquid species can lead to a weak modulation of
contact line friction upon changes in concentrations, whereas viscosity can be
affected much more substantially. It is not unreasonable to also envision the
opposite scenario, whereby the addition of molecular species with high
affinity with the substrate would increase contact line friction without
producing any dramatic change of bulk rheology. Such behaviors could be
desirable in applications where the time scale of wetting is shorter than the
possibly long relaxation times of polymers or diffusion times of colloidal
particles, such as droplet impact or rapid capillary rise.

##  VI Conclusions

We have presented the results of Molecular Dynamics simulations aimed to
investigate the correlation between contact line friction and viscosity of
water-glycerol droplets spreading over flat silica-like surfaces. The shear
viscosity coefficient has been computed from equilibrium and non-equilibrium
simulations and has been found to reflect the rheology of real water-glycerol
solutions. The contact line friction coefficient has been obtained directly by
fitting linearized Molecular Kinetic Theory to the contact angle and the
contact line speed measured from simulations. The scaling between contact line
friction and viscosity has been determined to be sub-linear, contrarily to the
predictions of single-specie MKT. A correction inspired by multicomponent MKT
has been proposed and linear scaling between interfacial friction and
viscosity has been recovered by accounting for the local change in mass
fraction and the difference of excluded volume between water and glycerol.

Even if the effect of molecular depletion can be noticeable on flat
substrates, simulating wetting dynamics on idealized surfaces constitutes a
significant simplification. In the case of aqueous glycerol droplets,
topographical roughness is known to affect the equilibrium contact angle on
silica surfaces [53, 54]. Therefore, simulating droplet spreading on more
realistic surfaces is a natural continuation of this work. Furthermore, the
assumptions introduced in section V.2 to correct the scaling between shear and
friction coefficients may be relaxed to address cases where to a change in
viscosity corresponds a substantial change in molecular kinetic rates, such as
when the wetting fluid or the surface are cooled down or heated up.

The most significant conclusion of this work is that molecular-scale effects,
which cannot be easily captured in a continuous description of liquid-solid
interfaces, contribute substantially to contact line friction. Even the rather
simple case of an aqueous solution wetting a flat surface shows deviation from
the postulated linear scaling. In essence, molecular heterogeneity emerges at
the interface scale and hence needs to be accounted for in order to correctly
characterize wetting phenomena.

## Supplementary information

See supplementary information at [URL will be inserted by publisher] for
additional details on uncertainty quantification, the measurement of molecular
interfacial observables and the calculation of shear viscosity.

## Data availability statement

All input of molecular simulations are available on Zenodo [55, 56, 57, 58,
59, 60]. Part of the output of is also available; very heavy files such as
fully-detailed trajectory have been excluded. Zenodo datasets also contain
simple Python scripts to open and visualize the output of molecular
simulations. More advanced analysis scripts can be provided under reasonable
request. The Gromacs fork for the density binning code is available on the
Github repository curated by Dr. Petter Johansson [61].

## Acknowledgements

We wish to thank Prof. Shervin Bagheri at KTH, Prof. Gustav Amberg at
Södertörns högskola and Dr. Petter Johansson at Université de Pau for the
useful feedback. Numerical simulations were performed on resources provided by
the Swedish National Infrastructure for Computing (NAISS snic2022-1-15) at
PDC, Stockholm. We acknowledge funding from Swedish Research Council
(INTERFACE Centre grant No. VR-2016–06119 and grant No. VR-2014–5680).

## Author contributions

Michele Pellegrino: conceptualization (equal), project administration (equal),
data curation, formal analysis, investigation, methodology, visualization,
writing - original draft, writing - review & editing (equal). Berk Hess:
conceptualization (equal), project administration (equal), supervision,
writing - review & editing (equal).

##  Appendix A Molecular Dynamics simulation

The parametrization of molecular simulations is analogous to the one of [62]
and [45]. We report the information needed to reproduce the simulations.
Configuration files can be found on Zenodo [55, 56, 57, 58, 59, 60].

Oxygen mass |  $m_{O}$ |  $9.95140$ | u
---|---|---|---
L-J well depth (silicon) |  $\varepsilon_{Si}$ |  $0.2$ | kJ mol-1
L-J char. distance (silicon) |  $\sigma_{Si}$ |  $0.45$ | nm
L-J well depth (oxygen) |  $\varepsilon_{O}$ |  $0.65019$ | kJ mol-1
L-J char. distance (oxygen) |  $\sigma_{O}$ |  $0.316557$ | nm
Si-O bond distance |  $d_{so}$ |  $0.151$ | nm
Hexagonal lattice spacing |  $d_{hex}$ |  $0.45$ | nm
Restraint force constant |  $\kappa_{O}$ |  $10^{5}$ | kJ mol-1 nm-2
Partial charge (oxygen) |  $q$ |  $-0.72$ | $|\mbox{e}^{-}|$
Table 5: Parametrization of the force field of silica quadrupoles from [62].

The force field of silica quadrupoles can be decomposed into a short-ranged
non-bonded van der Waals term, a long-range attractive electrostatic term and
a restrain term. Short-range non-bonded interactions are computed via the
Lennard-Jones potential
$U_{LJ}(r)=4\varepsilon[(\sigma/r)^{12}-(\sigma/r)^{6}]$ and long range
electrostatic via the Coulomb potential
$U_{C}(r)=(q_{i}q_{j})/(4\pi\epsilon_{0}r)$. Wall atoms are restrained to
fixed positions in space with a spring potential of the kind:
$U_{R}(r)=k(r-r_{ref})^{2}/2$. Restraining surface atoms instead of freezing
them in place allows momentum exchange between liquid and solid molecules,
which is necessary in order to avoid very large and unphysical slip lengths
[63]. Covalent bonds between oxygen and silicon atoms are treated as rigid
constraints. Silicon atoms are treated as virtual sites without mass. It has
been shown that the wetting behaviour on this type of surface is qualitatively
different from the one of simple Lennard-Jones crystals and it is better
suited to explore wetting dynamics on hydrophilic surfaces [13]. Table 5
reports the force field parameters for silica.

Lennard-Jones parameters for the interactions between silica and water or
glycerol are generated via the geometric combination rule. Short-range non-
bonded interactions are computed using a Verlet list and 1 nm cutoff distance
in real space. Electrostatic interactions are computed using the Particle Mesh
Ewald method (PME), with 1 nm cutoff distance for the calculation in real
space and 0.15 nm grid spacing for the calculation in the reciprocal space.

Integration over time is performed using the leapfrog algorithm with a time
step of 2 fs. Systems are thermalized using the Bussi-Donadio-Parrinello
thermostat (v-rescale in Gromacs), with a 1 ps coupling time. Ideal purely-
repulsive Lennard-Jones walls are placed beyond the silica surfaces at the
location of periodic boundary conditions along the vertical direction ensuring
no interaction between periodic images along $z$.

##  Appendix B Interface extraction

The measurement of the contact line speed and contact angle requires samples
of the position of the liquid-vapor interface over time. Interface extraction
for configurations of type DROP I and DROP II follows the procedure already
presented in [62] and [45].

Two-dimensional ($xz$) density maps for the liquid-vapor systems are computed
by binning the atomic positions on-the-fly, that is concurrently with the
running simulation, on a structured grid with spacing $\Delta x=\Delta z\simeq
0.2$ nm. The position of the interface along the vertical direction is in turn
obtained by linearly interpolating the center of the two cells where the
density is closer to half the bulk liquid density.

The contact angle is obtained in two different ways depending on the
configuration and the goal. When running configurations of type DROP I the
contact angle is identified via the tangent to the best least-squares circle
fit of the interface shape obtained after a sufficiently long time. On the
other hand, when simulating configurations of type DROP II we are interested
in the non-equilibrium transient; in this case the contact angle is obtained
by linearly interpolating the first 3 interface points, starting from 0.4 nm
above the reference position of silicon atoms.

Figure 6: Relaxation of the contact line position and the dynamic contact
angle over time. Gray lines and dots represent the results extracted directly
from molecular simulations without post-processing (fitting for the contact
line position and replica averaging for the contact angle).

The speed of the contact line is obtained from the displacement of the contact
line from its position at the beginning of the simulation. Due to thermal
fluctuations, the position of the contact line is not differentiable over
time; hence it is first fitted with a rational polynomial:

$x_{cl}(t)=\frac{a_{3}t^{3}+a_{2}t^{2}+a_{1}t+a_{0}}{a_{-1}t+1}\;,$ (16)

and then differentiated analytically [30]. Constraints on the coefficients
$a_{j}$ are imposed to ensure that the velocity is always positive and that
the acceleration is always negative, as expected in an overdamped spreading
regime with no inertial oscillations. On the other hand, no additional
filtering is performed on the signal of the contact angle over time except
replica averaging. Figure 6 shows the relaxation curves of the contact line
position and the dynamic contact angle over time for the configurations of
type DROP II.

##  Appendix C Measurement of shear viscosity

In this appendix we briefly describe how viscosity is computed from molecular
simulations. A thorough explanation of equilibrium and non-equilibrium methods
for the calculation of transport coefficients is outside the scope of this
article. Nonetheless, this appendix, the additional details reported in the
supplementary information and the data in [55, 56, 57, 58, 59, 60] ensure the
reproducibility of our calculations.

We first illustrate how viscosity is obtained using the equilibrium approach.
Systems of type BOX I are first equilibrated at constant temperature (300 K)
and pressure (1 bar) for 4 ns, after which long simulations ($>$100 ns) are
performed at constant temperature and volume. The ensemble average is
performed over 60 replicas, being either parallel runs with different initial
conditions sampled from the equilibration run or different slices of the same
molecular dynamics trajectory. The viscosity coefficient is obtained by linear
regression of the ensemble average against time. Given the diffusive nature of
the observable, the regression is weighted on the estimated standard error at
a given time [64]. The regression is performed within a time window determined
on one hand by convergence to linear behavior (visually inspected) and on the
other by a 10% threshold on the standard error.

The viscosity coefficient obtained using the non-equilibrium approach is
computed by fitting the velocity field binned from molecular trajectory to
expression 7. As briefly explained in section III.2, a drawback of this
approach is that a few values of the intensity of the external field $\xi$
need to be screened when comparing the resulting viscosity with the one
obtained with the equilibrium approach. In principle, a very small value of
$\xi$ may be desirable in order to remain consistent with the ‘linear
response’ assumption. However, the signal-to-noise ratio on the flow field is
better the larger the intensity of the perturbation. We considered 3
logarithmically spaced values for $\xi=\\{10^{-2},10^{-3},10^{-4}\\}$ nm/ps2.
Configurations of type BOX II equilibrated for 4 ns at constant temperature
and pressure. The non-equilibrium simulations are then performed at constant
volume and last at least 4 ns or until less than 10% relative error (computed
via bootstrapping) is achieved on the estimated viscosity.

##  Appendix D Interfacial friction

In this last appendix we report additional information on the coefficients
reported in table 4, used to evaluate equation 12. The lifetime of hydrogen
bonds between liquid and solid molecules is computed using the Gromacs command
gmx hbond -life on 3 ns long trajectories sampled every 1 ps. The per-molecule
volume of water and glycerol can be computed indirectly by running gmx
freevolume and supplying to the -select flag the group name of water or
glycerol molecules. The per-molecule volume is then computed knowing the total
number of molecules of one species and the difference in excluded volume
between the other species and the total system. The local mass fraction of
glycerol is computed in the layering region of MENISCUS systems within 0.8 nm
from the silica wall and sufficiently far away from the liquid-vapor interface
(as shown in figure 4) using on-the-fly binned density maps. The effective
local viscosity is computed by using equation 8 with $\alpha^{*}_{g}$ instead
of $\alpha_{g}$.

## References

  * Weinstein and Ruschak [2004] S. J. Weinstein and K. J. Ruschak, Annual Review of Fluid Mechanics 36, 29 (2004), https://doi.org/10.1146/annurev.fluid.36.050802.122049 .
  * de Ruiter  _et al._ [2017] R. de Ruiter, P. Colinet, P. Brunet, J. H. Snoeijer, and H. Gelderblom, Phys. Rev. Fluids 2, 043602 (2017).
  * Bureš and Sato [2022] L. Bureš and Y. Sato, Journal of Fluid Mechanics 933, A54 (2022).
  * Solga  _et al._ [2007] A. Solga, Z. Cerman, B. F. Striffler, M. Spaeth, and W. Barthlott, Bioinspiration & Biomimetics 2, S126 (2007).
  * Mouterde  _et al._ [2017] T. Mouterde, G. Lehoucq, S. Xavier, A. Checco, C. T. Black, A. Rahman, T. Midavaine, C. Clanet, and D. Quéré, Nature Materials 16, 658 (2017).
  * Voinov [1976] O. V. Voinov, Fluid Dynamics 11, 714 (1976).
  * Cox [1986] R. G. Cox, Journal of Fluid Mechanics 168, 169 (1986).
  * Petrov and Petrov [1992] P. Petrov and I. Petrov, Langmuir 8, 1762 (1992).
  * Shikhmurzaev [1997] Y. D. Shikhmurzaev, Journal of Fluid Mechanics 334, 211 (1997).
  * Blake [2006] T. D. Blake, Journal of Colloid and Interface Science 299, 1 (2006).
  * Do-Quang  _et al._ [2015] M. Do-Quang, J. Shiomi, and G. Amberg, EPL (Europhysics Letters) 110, 46002 (2015).
  * Carlson  _et al._ [2012a] A. Carlson, G. Bellani, and G. Amberg, Phys. Rev. E 85, 045302 (2012a).
  * Johansson  _et al._ [2015] P. Johansson, A. Carlson, and B. Hess, Journal of Fluid Mechanics 781, 695–711 (2015).
  * Johansson and Hess [2018] P. Johansson and B. Hess, Phys. Rev. Fluids 3, 074201 (2018).
  * Blake and Haynes [1969] T. Blake and J. Haynes, Journal of Colloid and Interface Science 30, 421 (1969).
  * Blake and De Coninck [2002] T. Blake and J. De Coninck, Advances in Colloid and Interface Science 96, 21 (2002), a Collection of Papers in Honour of Nikolay Churaev on the Occasion of his 80th Birthday.
  * Blake and Shikhmurzaev [2002] T. Blake and Y. Shikhmurzaev, Journal of Colloid and Interface Science 253, 196 (2002).
  * Ramiasa  _et al._ [2013] M. Ramiasa, J. Ralston, R. Fetzer, R. Sedev, D. M. Fopp-Spori, C. Morhard, C. Pacholski, and J. P. Spatz, Journal of the American Chemical Society 135, 7159 (2013).
  * Andrukh  _et al._ [2014] T. Andrukh, D. Monaenkova, B. Rubin, W. Lee, and K. G. Kornev, Soft Matter 10, 609 (2014).
  * Mohammad Karim  _et al._ [2016] A. Mohammad Karim, S. H. Davis, and H. P. Kavehpour, Langmuir 32, 10153 (2016), pMID: 27643428, https://doi.org/10.1021/acs.langmuir.6b00747 .
  * Barrio-Zhang  _et al._ [2020] H. Barrio-Zhang, É. Ruiz-Gutiérrez, S. Armstrong, G. McHale, G. G. Wells, and R. Ledesma-Aguilar, Langmuir 36, 15094 (2020).
  * Duvivier  _et al._ [2011] D. Duvivier, D. Seveno, R. Rioboo, T. D. Blake, and J. De Coninck, Langmuir 27, 13015 (2011).
  * Duvivier  _et al._ [2013] D. Duvivier, T. D. Blake, and J. De Coninck, Langmuir 29, 10132 (2013), pMID: 23844877, https://doi.org/10.1021/la4017917 .
  * Seveno  _et al._ [2009] D. Seveno, A. Vaillant, R. Rioboo, H. Adão, J. Conti, and J. De Coninck, Langmuir 25, 13034 (2009).
  * Chen  _et al._ [2014] L. Chen, J. Yu, and H. Wang, ACS Nano 8, 11493 (2014), pMID: 25337962, https://doi.org/10.1021/nn5046486 .
  * Deng  _et al._ [2016] Y. Deng, L. Chen, Q. Liu, J. Yu, and H. Wang, The Journal of Physical Chemistry Letters 7, 1763 (2016), pMID: 27115464, https://doi.org/10.1021/acs.jpclett.6b00620 .
  * Carlson  _et al._ [2012b] A. Carlson, G. Bellani, and G. Amberg, EPL (Europhysics Letters) 97, 44004 (2012b).
  * Li  _et al._ [2023] X. Li, F. Bodziony, M. Yin, H. Marschall, R. Berger, and H.-J. Butt, Nature Communications 14, 4571 (2023).
  * Papadopoulou  _et al._ [2019] E. Papadopoulou, C. M. Megaridis, J. H. Walther, and P. Koumoutsakos, ACS Nano 13, 5465 (2019).
  * Fernández-Toledano  _et al._ [2020] J. C. Fernández-Toledano, T. D. Blake, J. De Coninck, and M. Kanduč, Phys. Rev. Fluids 5, 104004 (2020).
  * Ozsipahi  _et al._ [2022] M. Ozsipahi, Y. Akkus, and A. Beskok, International Communications in Heat and Mass Transfer 136, 106166 (2022).
  * Giri  _et al._ [2022] A. K. Giri, P. Malgaretti, D. Peschka, and M. Sega, Phys. Rev. Fluids 7, L102001 (2022).
  * Takamura  _et al._ [2012] K. Takamura, H. Fischer, and N. R. Morrow, Journal of Petroleum Science and Engineering 98-99, 50 (2012).
  * Yada  _et al._ [2023] S. Yada, K. Bazesefidpar, O. Tammisola, G. Amberg, and S. Bagheri, Phys. Rev. Fluids 8, 043302 (2023).
  * Liang  _et al._ [2010] Z.-P. Liang, X.-D. Wang, Y.-Y. Duan, Q. Min, C. Wang, and D.-J. Lee, Langmuir 26, 14594 (2010).
  * Collins [2004] F. C. Collins, The Journal of Chemical Physics 26, 398 (2004), https://pubs.aip.org/aip/jcp/article-pdf/26/2/398/11083630/398_1_online.pdf .
  * Abraham  _et al._ [2015] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1-2, 19 (2015).
  * Jorgensen  _et al._ [1996] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, Journal of the American Chemical Society 118, 11225 (1996).
  * Jahn  _et al._ [2014] D. A. Jahn, F. O. Akinkunmi, and N. Giovambattista, The Journal of Physical Chemistry B 118, 11284 (2014).
  * Hess [2002] B. Hess, The Journal of Chemical Physics 116, 209 (2002), https://aip.scitation.org/doi/pdf/10.1063/1.1421362 .
  * Walther  _et al._ [2013] J. H. Walther, K. Ritos, E. R. Cruz-Chu, C. M. Megaridis, and P. Koumoutsakos, Nano Letters 13, 1910 (2013).
  * Volk and Kähler [2018] A. Volk and C. J. Kähler, Experiments in Fluids 59, 75 (2018).
  * Segur and Oberstar [1951] J. B. Segur and H. E. Oberstar, Industrial & Engineering Chemistry 43, 2117 (1951), https://doi.org/10.1021/ie50501a040 .
  * Cheng [2008] N.-S. Cheng, Industrial and Engineering Chemistry Research 47, 3285 (2008), https://doi.org/10.1021/ie071349z .
  * Pellegrino and Hess [2022] M. Pellegrino and B. Hess, Physics of Fluids 34, 102010 (2022), https://doi.org/10.1063/5.0121144 .
  * Cheng  _et al._ [2001] L. Cheng, P. Fenter, K. L. Nagy, M. L. Schlegel, and N. C. Sturchio, Phys. Rev. Lett. 87, 156103 (2001).
  * Bitsanis  _et al._ [1987] I. Bitsanis, J. J. Magda, M. Tirrell, and H. T. Davis, The Journal of Chemical Physics 87, 1733 (1987), https://pubs.aip.org/aip/jcp/article-pdf/87/3/1733/13395951/1733_1_online.pdf .
  * Dashnau  _et al._ [2006] J. L. Dashnau, N. V. Nucci, K. A. Sharp, and J. M. Vanderkooi, The Journal of Physical Chemistry B 110, 13670 (2006), pMID: 16821896, https://doi.org/10.1021/jp0618680 .
  * Hoang and Galliero [2012] H. Hoang and G. Galliero, Phys Rev E Stat Nonlin Soft Matter Phys 86, 021202 (2012).
  * Joanny and de Gennes [1984] J. F. Joanny and P. G. de Gennes, The Journal of Chemical Physics 81, 552 (1984), https://pubs.aip.org/aip/jcp/article-pdf/81/1/552/9722371/552_1_online.pdf .
  * Perrin  _et al._ [2016] H. Perrin, R. Lhermerout, K. Davitt, E. Rolley, and B. Andreotti, Phys. Rev. Lett. 116, 184502 (2016).
  * Lu  _et al._ [2016] G. Lu, X.-D. Wang, and Y.-Y. Duan, Advances in Colloid and Interface Science 236, 43 (2016).
  * Hitchcock  _et al._ [1981] S. J. Hitchcock, N. T. Carroll, and M. G. Nicholas, Journal of Materials Science 16, 714 (1981).
  * Moldovan  _et al._ [2014] A. Moldovan, M. Bota, D. Dorobantu, I. Boerasu, D. Bojin, D. Buzatu, and M. Enachescu, Journal of Adhesion Science and Technology 28, 1277 (2014), https://doi.org/10.1080/01694243.2014.900907 .
  * Pellegrino [2023a] M. Pellegrino, Surface tension coefficient of acqueous glycerol from Molecular Dynamics simulations, https://doi.org/10.5281/zenodo.8024992 (2023a).
  * Pellegrino [2023b] M. Pellegrino, Shear viscosity coefficient of acqueous glycerol from equilibrium Molecular Dynamics simulations, https://doi.org/10.5281/zenodo.7576495 (2023b).
  * Pellegrino [2023c] M. Pellegrino, Shear viscosity coefficient of acqueous glycerol from non-equilibrium Molecular Dynamics simulations, https://doi.org/10.5281/zenodo.7756756 (2023c).
  * Pellegrino [2023d] M. Pellegrino, Equilibrium contact angle of acqueous glycerol on a silica-like surface from Molecular Dynamics, https://doi.org/10.5281/zenodo.8027883 (2023d).
  * Pellegrino [2023e] M. Pellegrino, Molecular Dynamics simulations of of acqueous glycerol spreading on a silica-like surface, https://doi.org/10.5281/zenodo.8047708 (2023e).
  * Pellegrino [2023f] M. Pellegrino, Equilibrium Molecular Dynamics of aqueous glycerol menisca, https://doi.org/10.5281/zenodo.8077915 (2023f).
  * Johansson [2022] P. Johansson, Fork of gromacs used to perform flow variables binning, https://github.com/pjohansson/gromacs-flow-field.git (2022), [GitHub repository].
  * Lācis  _et al._ [2022] U. Lācis, M. Pellegrino, J. Sundin, G. Amberg, S. Zaleski, B. Hess, and S. Bagheri, Journal of Fluid Mechanics 940, A10 (2022).
  * Sega  _et al._ [2013] M. Sega, M. Sbragaglia, L. Biferale, and S. Succi, Soft Matter 9, 8526 (2013).
  * Zhang  _et al._ [2015] Y. Zhang, A. Otani, and E. J. Maginn, Journal of Chemical Theory and Computation 11, 3537 (2015).