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

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


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
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 θ is allowed to deviate from its equilibrium value θ 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: µ f being referred to as the contact line friction coefficient and σ being the liquid-vapor interface tension; λ is the length of molecular jumps and κ 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 ∆g according to Arrhenius equation: The jump activation free energy can be split into one component that quantifies the adsorption (desorption) work on (respectively from) the solid surface ∆g s , and one encompassing the interactions between neighboring liquid molecules ∆g l [16].The latter is in turns related to the viscosity of the wetting fluid η according to Eyring theory: η ∝ exp{∆g l /(k B T )}/v, being v the volume of a fluid molecule [36].Hence, the equilibrium jump rate is decomposed as: Note that κ 0 s and κ 0 l 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: 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.

A. 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 α 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.
In order to fully characterize the bulk, interfacial and wetting properties of the waterglycerol-silica combination, several different configurations of bulk liquid, liquid-vapor interfaces and liquid-vapor-solid contact lines need to be simulated.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.

B. 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: 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]: 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: In contrast to the equilibrium approach, the viscosity could in principle depend on the amplitude of the external acceleration.Therefore, a few values of ξ need to be tried to assess the consistency of results.

A. Material parameters and transport coefficients
Table II 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 α 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).
The calculation of viscosity via equilibrium and non-equilibrium simulations is performed for each value of α g .The results are presented in figure 2a.Although it was not possible to Comparison between Molecular Dynamics results, experimental results by Segur and Oberstar [43] and the best fit of the empirical model by Cheng [44]; β is defined in equation 8.
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 α g and ξ.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]: 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.

B. Contact line friction
Contact line friction is estimated directly by fitting equation 1 to the post-processed results of molecular simulations, with fitting parameters {µ f , θ 0 } for each glycerol fraction.
The equilibrium contact angle can be obtained from the intercept of the linear fit between u cl and cos θ, or can be imposed to the value obtained separately from equilibrium runs (table II).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: being θ 0 the value of the equilibrium contact angle reported in table II.The notation ⃗ • indicates the list of values obtained from molecular simulations.The contact line friction parameter resulting from the minimization of L is found to be insensitive to the penalty term; the penalty coefficient ζ is set to 0.1.Further information on uncertainty quantification can be found in the supplementary information.
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 α g : this is a signature of contact line friction.Furthermore, the linear MKT model appears suitable to capture the correla- tion 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 III.
Note that no result is reported for α 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 III. Figure 3b  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 α 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: 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

A. 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 C, 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: being n the number of molecular species.The superscript • * 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 α * g < α g , owing to the results of section IV C. The dependence on viscosity is 'hidden' in the equilibrium jump rate κ 0 i ∝ κ 0 s,i /(v i η * ), 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.

B. Correction to linear scaling
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 C 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 λ w ≃ λ g ≃ 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: We now introduce a further assumption: η * = η(α * 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 α * 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 IV.Appendix D briefly explains how the coefficients are obtained from molecular simulations.
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 leadingorder correction for depletion and excluded volume effects is therefore enough to reconcile the results of Molecular Dynamics simulation with Molecular Kinetic Theory.

C. 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: 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 µ app f and µ f [8]: 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 θ(u cl ) and 14 into 13 for the apparent contact angle.After linearizing for small wetting speeds we are left with the following expression: 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.
D. 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 µ f ≫ η, where wetting is thus substratedominated [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.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 B 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.
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].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 V 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.
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: 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.The SASA (solvent accessible surface area) is computed from systems of type SLAB using the command gmx sasa; we refer to Gromacs documentation for a detailed explanation of the command functionality [2].
The calculation of surface tension may be inaccurate due to the long-range treatment of dispersion forces or due to limitations in the water model.We have quantified how sensitive the estimate of contact line friction (from equation 1 in the main article) is to the change of surface tension, assuming the spreading dynamics remains the same and fixing either (i) the equilibrium contact angle, or (ii) the difference between solid-vapour and solid-liquid surface energies (σ SV − σ SL ).We then estimated the error in the calculation of contact line friction using: where σ is either the experimental surface tension or (a) the one computed using TIP4P/2005    in order to correctly sample the velocity autocorrelation function in time.Since we are interested in the long-time behavior of the integral it is unnecessary to output energy every step and it is sufficient to set nstenergy=100.

FIG. 1 .
FIG. 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. Scheme

FIG. 2 .
FIG. 2. (a): Shear viscosity (in centiPoise) as a function of the non-equilibrium forcing amplitude ξ and the mass fraction of glycerol α 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):

FIG. 3 .
FIG. 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].
FIG. 4. (a): Semi-logarithmic plot of the density profiles of water and glycerol as a function of the vertical coordinate, for the case α 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.

FIG. 5 .
FIG. 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 α * g , see supplementary information.
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.

Figure 3
FIG. 3: Surface tension of aqueous glycerol solutions for SPC/E and TIP4P/2005 with OPLS-AA glycerol, compared to experimental results from [3].

and 1 nm
FIG. 4: Sensitivity of the scaling between contact line friction and viscosity due to the effect of real-space Lennard-Jones cutoff on surface tension.Left: case (i), right:case (ii); see equation 1.

FIG. 7 :
Figures 6a and 7a report the integral returned by gmx energy -eviscoi for each replica and its average, for α g = 0 and α g = 0.8.It can be noted how the Brownian behavior of the observable leads to more uncertain estimates for longer times.We compute the relative standard error for each fixed time frame across 60 replicas and impose an arbitrary threshold of 7.5%.

TABLE II .
9±0.1 5.59±0.036.14±0.1848.94±1.60Density, 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 ( †) is not reported for TIP4P/2005 since it does not depend on the water model.

TABLE III .
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.

TABLE IV .
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.

TABLE I :
Liquid-vapour surface tension for different glycerol concentrations and differentLennard-Jones cutoff radii, compared to PME Lennard-Jones and experimental results.
FIG. 2: Left: visualization of a system of type SLAB.Right: solvent accessible surface area of water and glycerol w.r.t. the interface with vapor.

TABLE II :
Results of the uncertainty quantification procedure (average ± standard deviation) for each of the model parameters, for different weight values.