Asymmetry of wetting and de-wetting on high-friction surfaces originates from the same molecular physics

The motion of fluid–fluid interfaces over solid substrates is a common occurrence in processes both related to the natural world and to industry.^1,2 First attempts to describe the dynamics of a contact line (c.l.) using the tools of continuous fluid mechanics encountered an obstacle represented by the fact that no-slip boundary conditions entail infinite stresses at the substrate, would the contact line move.³ Chasing this mathematical paradox, many different attempts to fully regularize the problem have emerged.^4–7 Among those, one of the most popular is the contact line friction model which stems from the Molecular Kinetic Theory (MKT) by Blake and Haynes:⁸ under the assumptions of negligible inertial effects and total dissipation of the fluid's kinetic energy due to friction localized in the contact line region, one obtains a direct mobility relation between contact line speed and microscopic dynamic contact angle: $Ucl(t)=f(θ(t))$⁠.

While the existence of a microscopic dynamic contact angle has now been widely accepted by the fluid dynamic and wetting communities, the existence and the characterization of the speed-angle relation still remains controversial and criticized.⁹ One of the key arguments against the employment of mobility models is the difficulty to reconcile results obtained under spontaneous wetting, where contact lines are driven by capillary action due to out-of-equilibrium initial conditions, and forced wetting, where contact lines are driven by mechanical action.¹⁰ Proving or disproving contact line models via physical experiments is inherently arduous due to the optical resolution required to resolve the interface curvature in the contact line region. Recently, studies involving optical microscopy and atomic force microscopy have managed to infer the value of microscopic dynamic contact angles.^11–13 However, finely controlling flow conditions and surface wettability in experiments still remains challenging. From this perspective, atomistic Molecular Dynamics (MD) simulations can be used in lieu of experiments to access nanoscopic length and time scales and obtain almost unrestrained control over fluid and surface features. Both spontaneous and forced wetting modes are easily reproducible in MD by either placing a droplet on a surface and changing the surface wettability so that the initial contact angle does not match its equilibrium value (spontaneous), or by confining it between two solid walls and shear them in opposite directions to reproduce a two-phase Couette flow (forced).

Despite MD gaining popularity in studying contact line problems, most research is focused on rather simple cases, such as Lennard–Jones fluids wetting Lennard–Jones crystals. Simulating more complex fluid/surface combinations usually demands more elaborate molecular models and more computational resources. Despite technical challenges, recent works have shown how non-equilibrium MD (NEMD) simulations of water droplets confined between silica-like walls can be effectively employed to characterize continuous simulation techniques “equipped” with dynamic wetting boundary conditions that do not rely on imposing unphysically large slip lengths.^14,15 Furthermore, upon directly inspecting the dynamic contact angle from NEMD simulations and correlating it with the contact line velocity, it was discovered that advancing and receding contact lines experience different contact line friction; this phenomenon has not been observed for Lennard–Jones substrates.¹⁶ In particular, the contact line friction on hydrophobic substrates is smaller for receding contact lines and larger for advancing ones, while hydrophilic substrates show the opposite behavior.¹⁵ Interpreting the motion of advancing c.l. as “wetting” and the one of receding c.l. as “de-wetting,” hydrophilic substrates are easier to wet rather than de-wet, and vice versa. This seemingly intuitive observation fueled our interest to further investigate the asymmetry between advancing and receding contact lines.

MKT predicts an odd functional relationship between contact line speed and uncompensated Young's stress, i.e., it predicts the same contact line speed for the same contact angle deviation, with the sign corresponding to whether the contact line is advancing or receding; hence, it cannot capture contact line friction asymmetry naturally. A recent model inspired by the physico-chemistry of water on high-friction surfaces predicts the contact line coefficient to change depending on the absolute value of the dynamic contact angle¹⁷ and, thus, represents a good candidate to automatically capture contact line asymmetry. While MKT considers only molecular motion in the wall-nearest liquid layer responsible for the motion of the contact line, this new model adds the contribution of transport from the liquid layer above (in the form of “rolling” motion). It is well known that hydrophilic polar surfaces modify the hydration structure of interfacial water;¹⁸ the importance of considering two fluid layers instead of one is corroborated by simple observations on the conformation of water molecules in the contact line region. As an additional effort toward showing the importance of hydration structure, we propose a method to systematically connect molecular conformation to molecular motion in order to extract kinematic information regarding the processes that advance or recede contact lines. The method involves analyzing molecular transitions with the toolbox of Markov state models.

The purposes of this work are to advance the physical understanding of contact line friction and to show how Molecular Dynamics can be employed to characterize the mobility of three-phase contact lines. These newfound insights can lead to a better modeling of the boundary conditions of continuous fluid mechanics. In order to ensure the reproducibility of our work and to share a set of MD benchmarks useful to further validate contact line models, the output of our MD simulations will be provided in openly accessible databases.

This paper is organized as follows: in Sec. II, we illustrate the mathematical model of contact line friction. In Sec. III, we illustrate the NEMD and the analysis techniques employed in this work. In Sec. IV, we illustrate and briefly discuss the main results of molecular simulations. In Sec. V, we present the study of molecular kinematics of contact lines. In Sec. VI, we further discuss the implications of our results, and finally, in Sec. VII, we draw the conclusions about the impact of the present work and how to proceed forward.

We consider the case of a two-dimensional droplet or meniscus. The wetted surface is atomistically flat and chemically homogeneous. Under such conditions, the equilibrium wetting configuration of a droplet or a meniscus can be uniquely identified by the Young–Dupré contact angle θ₀ between the liquid–vapor interface and the solid surface [Fig. 1(a)].

The contact line friction model assumes that when a contact line moves with steady velocity U_cl, the microscopic contact angle θ deviates from its equilibrium value, causing an uncompensated Young's stress, quantifiable as $FY=γ(cos θ0−cos θ)$⁠, with γ being the liquid–vapor surface tension. Molecular Kinetic Theory describes the motion of contact lines as a sequence of thermally activated molecular “jumps:” liquid molecules will preferably remain close to the minima of the potential energy surface of the solid substrate and rattle until they have enough thermal energy to escape, jump, and get adsorbed in another local minimum [Fig. 1(b)]. The mobility of contact lines is, therefore, determined by the “corrugation” of the surface energy landscape rather than the equilibrium contact angle itself.¹⁹ Furthermore, thermally activated molecular hopping is inherently dissipative. Hence, MKT predicts a direct mobility relation between U_cl and F_Y, which after linearizing for small angle deviations reads as

where μ_f is the contact line friction parameter. Although contact line friction depends ultimately on the details of the potential energy surface, it can be correlated with the equilibrium contact angle if the surface topography and chemical composition are kept fixed.^15,16

Equation (1) is symmetric in the sense that it predicts a moving contact line to experience the same friction for the same absolute speed (or equivalently the same deviation from the equilibrium contact angle), regardless of whether the contact line is advancing (i.e., it moves toward the liquid phase) or receding (i.e., it moves toward the vapor phase). The same consideration holds for the full non-linear mobility expression. MKT is based on the assumption that molecules adsorbed in the potential minima of the substrate are mobile enough to effectively carry the contact line motion. If the surface is very hydrophilic, however, it may happen that energy barriers in this adsorption/desorption process are too high and it becomes energetically preferable to advance a contact line by transporting molecules from the above liquid layers instead [Fig. 2(a)]. Johansson and Hess¹⁷ included this phenomenology in the contact line friction model by envisioning an angle-dependent contact line friction parameter,

The reasoning behind Eq. (2) is that molecules in the second layer will find the harder to “roll over” the stronger the interactions with neighboring water molecules and the smaller the dynamic contact angle. Parameter a captures the effect of the neighboring liquid molecules; $ε(θ)$ is an angle-dependent energy barrier, which scales with distance needed for a molecule in the second fluid layer to roll over an adsorbed molecule in the first layer. Such a distance will depend on the inclination angle between the two layers, that is, the dynamic contact angle. Using simple geometric arguments, the barrier has been previously estimated by Johansson and Hess¹⁷ as $ε(θ)=(0.5 sin θ+cos θ)2$⁠. This contact line friction model (which we will refer to as two-layer model) not only incorporates a richer phenomenology compared to MKT but also predicts situations where friction is asymmetric, that is, the same contact line speed causes different deviation of the contact angle depending on the wetting direction [Fig. 2(b)].

When discussing results, the following non-dimensional quantities for the contact line friction and the contact line speed will be used,

with μ being the viscosity of the fluid. Since μ_f and μ have the same physical dimensions, one can interpret $μf*$ in terms of “how strong” is solid–liquid friction compared to the friction between liquid layers. $Cacl$ can be regarded as the contact line capillary number. Superscripts will be dropped for the sake of readability. Advancing and receding contact lines will be implicitly distinguished by the sign of the contact line capillary number (⁠$Cacl>0$⁠: advancing, $Cacl<0$⁠: receding).

The molecular systems are simulated using GROMACS. The liquid quasi-2D droplets [Fig. 3(a)] and menisci [Fig. 3(b)] consists, respectively, of 300 305 and 172 933 SPC/E water molecules.²⁰ The substrates are composed of silica quadrupole molecules arranged in an hexagonal lattice [Fig. 4(a)]. The substrate wettability can be tuned by changing the value of the partial charge q on the oxygen atoms of silica [Fig. 4(b)]. Albeit the lattice of silica quadrupoles is not physically realistic, the liquid-substrate combination is the simplest one reproducing the hydrophilic interaction that would characterize an high-friction surface [Fig. 4(c)].

Non-equilibrium MD simulations are driven either by the initial condition being out of equilibrium or by a steady external action. In droplet spreading simulations, a thin cylinder of SPC/E molecules is initialized so that the distance between the lower liquid–vapor interface and the upper oxygen atoms of the silica molecules is $∼0.7$ nm. The droplet will spontaneously wet the surface until an equilibrium state characterized by the contact angle $θ0(q)$ is reached. Then, by increasing q, the droplet will spread even more over the surface (which is made more hydrophilic), while upon lowering q, the droplet will retract and de-wet the surface (which is made more hydrophobic). Changing the wettability of the surface thus allows us to simulate both advancing and receding contact lines.

Shear droplet simulations, on the other hand, are driven by constantly moving the silica walls. The confined droplet is first initialized as a rectangular prism with height-to-width ratio $∼0.75$ and let relax until it reaches its equilibrium meniscus shape. Then, the silica walls are moved in opposite directions with velocity U_w. If the wall velocity is below a critical velocity $Ucr(θ0)$⁠, then the motion of the contact line is steady, i.e., the contact line velocity is zero if viewed from an absolute frame of reference or $−Uw$ if viewed from the frame of reference of the moving wall. Thus, the speed at which the contact line moves can be controlled. For $Uw≥Ucr$ a critical transition occurs, which depending on the surface's wettability manifests either as a liquid film deposition at the receding contact line or as vapor entrainment at the advancing contact line.

The speed of the contact line will always be defined with respect to the inertial frame of reference moving alongside the corresponding substrate. Viscosity and surface tension for SCP/E water are the same as the ones in Lācis et al.,^14,15 that is $γ=5.78×10−2$ Pa $·$ m and $μ=8.77×10−4$ Pa $·$ s. The critical capillary number for shear droplet configurations and the equilibrium contact angle as function of the partial charge q are reported in Table I. Additional information on the parametrization of the force field can be found in the supplementary material.

Given the size of the molecular system, analyzing fully detailed and finely sampled atomistic trajectories over scales of tens of nanoseconds would be prohibitive. The measurement of the contact line speed and contact angle does not require the knowledge the positions and velocities of all atoms at every time, but only of the position of the vapor–liquid interface; this one can be extracted as a contour line of the liquid density map $ρ(x,z)$⁠. The density map is obtained by binning the atomic positions on-the-fly, that is concurrently with the running simulation, on a structured grid with spacing $Δx=Δz≃0.2$ nm.

The extraction of the interface is performed using a local half-density criterion: for each coordinate z_i corresponding to the center of a bin the local density profile $ρj(x)$ is extracted; the local bulk density $ρj,bulk$ is defined as average of the density values in a region $Δxbulk=10$ nm, centered on the center of mass of the density slice [Fig. 5(a)]. The interface point then is obtained by linearly interpolating the center of the two cells where the density is closer to $1/2·ρj,bulk$⁠. This approach produces a smooth interface curve despite the weak density layering occurring close to solid walls [Fig. 5(b)]. Moreover, it is not required to measure the interface location for each z_j, but only for a few bins close to the solid substrate. Based on the analysis in Ref. 15, the contact line location is identified with the interface point in the bin at $zj=0.9$ nm (respectively $zj=29.466$ nm for the upper wall in the shear droplet setup), which is about 0.3 nm above the vertical coordinate of silicon atoms.

In an equilibrium configuration, one may identify the contact angle via the tangent to the best least squares circle fit of the interface shape, following Young–Laplace law. The same approach obviously cannot be used when out of equilibrium; in such cases, which are the most interesting for this work, a linear interpolation of the first three interface points is used instead. Among all the interpolation techniques tried, this one provided more consistent results with the results in Ref. 15.

Contact line velocity is not constant for spreading droplets. Obtaining the velocity from the signal of the contact line position over time via finite difference would yield extremely noisy results; hence, a rational polynomial fit is performed on the position and the speed is obtained via analytical differentiation, as in Ref. 16. A numerical scanning on the number of polynomial coefficients is performed in order to single out the functional forms that produces the best fit for advancing and receding contact lines,

The number of coefficients of the rational polynomial function is different depending if the contact line is advancing or receding, but it is kept the same for all equilibrium contact angles. It is important to remark that equations (4) carry no physical meaning and are employed only for filtering purposes.

In order to study the near-contact-line conformation and mobility or water molecules, we analyze detailed molecular trajectories sampled every 1ps in a time window of 4 ns. Since the full molecular system would be far too large to process, we select atomic positions from a control volume of size $lx≃2.82$ nm $× ly≃4.68$ nm $× lz≃1.44$ nm centered at the moving contact line. In shear droplet simulations, the contact line position is stationary in the frame of reference of the simulation box, making the selection easier. We focus on the case of $θ0=69°$ and $θ0=38°$⁠, at equilibrium.

Molecular coordinates are projected on the (x, y) plane and a clustering procedure using the DBSCAN algorithm is conducted to discriminate between the water molecules in the liquid wedge and the ones adsorbed in the substrate ahead of the contact line. Then, the contact line curve $xcl(y)$ is approximated by the subset of the concave hull (alpha-shape) facing the solid–vapor interface [Fig. 6(a)]. The contact line molecules are defined as the ones having the x-coordinate of oxygen atoms within $dH2O=0.3166$ nm from the contact line, $|xO,i−xcl(yO,i)|<dH2O$⁠. It is worth mentioning that this procedure assumes the molecules closer to the contact line to be also closer to the solid wall; hence, it can be employed only if $θ<90°$⁠.

The conformation of water molecules is defined via two degrees of freedom: z_O, the unsigned distance along z between the coordinate of the oxygen atom and the coordinate of substrate's silicon atoms, and $φ$⁠, the polar angle of the molecule's dipole with respect to the unit vector along z [Fig. 6(b)]. The 2D histogram of the distribution of d.o.f. $ω(zO,φ)$ is computed using a simple binning procedure. The histogram is then transformed into an energy landscape via the Boltzmann formula: $E=−kBT log [ω(zO,φ)]$⁠. The wells in the profile of E represent the preferable conformations for water molecules.

Figure 7 displays the measurements of dynamic contact angle and contact line speed for spreading and shear droplet simulations. Despite the measurements from spreading droplet simulations showing significant variability, quantitative agreement between forced and spontaneous wetting is evident for $θ0=95°$ and $θ0=69°$⁠. In the case of $θ0=127°$ and $θ0=38°$⁠, both wetting modes provide consistent results, albeit the range of contact line velocity is different; as a matter of fact, it is easier to collect measurements from droplet spreading simulations if the contact line is advancing when the substrate is hydrophilic (vice versa if hydrophobic). The simulations, therefore, do not show any significant difference between forced and spontaneous wetting configuration in terms of contact angle—contact line speed relation. This is consistent with recent state-of-the-art molecular dynamics investigations¹⁶ and corroborates the assumption that forms the bedrock of the contact line friction model, according to which the contact line mobility relation should be invariant on the wetting mode at the molecular scale.

Once verified that both spontaneous and forced wetting modes produce consistent results, the next step is to fit the model described by Eq. (2). In order to reduce the arbitrariness of the functional fit, the parameters a and μ_f are fitted using separate criteria. In particular, since the model predicts a to be dependent on the liquid composition but not on the surface's wettability, only one value will be fitted encompassing all the MD results. On the other hand, the contact line friction depends on the strength of the electrostatic interactions between water and silica, which in turns unequivocally expresses the equilibrium contact angle; hence, a different value of contact line friction will be fitted depending on the surface wettability $μf(θ0)$⁠. The obtained values of a and $μf(θ0)$ are shown in Table II and show good agreement with the ones previously obtained in Ref. 17. The uncertainty range is estimated by performing 25 cross-validation steps where 20% of the dataset is removed upon estimating the parameters. Additional details regarding uncertainty quantification are provided in the supplementary material.

The results in this section show the capability of the two-layer model to capture all four possible combinations of moving contact lines: forced and spontaneous, advancing and receding. In Secs. IV B and V, we portray the molecular physics supporting the model.

Figure 8(a) shows the contour lines of the energy landscape $E(zO, cos φ)$ obtained after the binning procedure, for $θ0=69°$⁠. Two wells can be spotted, the deeper having minimum at $(zO≃0.2 nm,cos φ≃1)$ and the other having minimum for $(zO≃0.43 nm,cos φ≃−0.63)$⁠. While weak horizontal layering of water molecules has already been observed in Ref. 15, it can be noticed that each water layer shows a preferential orientation of the electrostatic dipole. Water molecules in the wall-nearest layer have the dipole pointing upward; their conformation thus corresponds to what has been labeled “adsorption” in Fig. 4(c). Water molecules in the layer above have their dipole mainly pointing downward, although the distribution of $cos φ$ is more outspread; their conformation corresponds to what has been labeled “hydrogen bonding.” The landscape for $θ0=38°$ [Fig. 8(b)] appears qualitatively similar.

The profile of $E(zO, cos φ)$ suggests that molecules in different layers of the fluid interact with the substrate in qualitatively different ways [Fig. 8(b)]. It is, therefore, possible to speculate that molecules adsorbed in layer 1 would find it more difficult to escape the energy well compared to the ones in layer 2 due to being caged into the silica lattice. From $E(zO, cos φ)$⁠, it is even possible to compute the energy barrier between layers 1 and 2 by taking the minimum free energy pathway. Using a simple Monte Carlo procedure, the energy barrier is roughly estimated as $ΔE=1.31 ± 0.014kBT$⁠; this value is not much outside the uncertainty range of parameter a, which, according to the discussion in Sec. II, scales the energy barrier associated with rolling motion. It must be pointed out that this barrier is modulated by the collective conformations of water molecules in the contact line region. Therefore, the value of the barrier computed from E may overestimate the actual energy barrier for a molecular displacement event of a moving contact line and, thus, should be considered only a zeroth order approximation. However, this approach to estimate displacement barriers is not dissimilar to the one taken by MKT, where single molecular jumps are considered instead of collective processes.

In order to thoroughly validate the two-layer model, the kinematics of near-contact-line molecules needs to be studied. Reasoning in the framework of MKT, the speed of a contact line is dictated by the frequency of molecular displacement events and their distance. For the situation depicted in this work, the multicomponent molecular-kinetic approach by Liang et al.²¹ can prove instrumental in mapping kinematic information to different conformations of water molecules. Considering the contribution of rolling and “jumping” molecules to the displacement of the contact line alike the one of two different species, we can rewrite the speed of the contact line as

with κ_i and $κi0$ being the non-equilibrium and the equilibrium displacement rates of species i, δ_i its displacement length and

the exponential tilting of the potential energy landscape occurring when then contact line is out of equilibrium, with w_i being the specific work of adhesion and x_i the local mass fraction of species i. Ultimately, comparing the exponential tilting for each species should provide quantitative information about their relative effect to contact line mobility. However, a direct computation of expression 6 is unfeasible due to the difficulty in characterizing the work of adhesion. Conversely, assuming both κ_i and $κi0$ are quantifiable, one can obtain the tilting as $expti=κi/κi0$⁠.

The method that we employ to extract the displacement rates from molecular simulation involves modeling the transition between different water conformations as a quasi-Markov process and compute the corresponding transition matrix. Let us label state 1 as the one of molecules adsorbed in the silica lattice, 2 as the one of molecule in the second layer that are hydrogen bonded and 3 the one of molecules that are away from the contact line, either in the bulk or at interfaces [Fig. 9(a)]. A rolling event would correspond to the transition from 2 to 1. In order to define a jump event we need to incorporate memory effects and define an additional state 1^* which is employed only to describe adsorption to a different silica lattice cavity. Therefore, the transition from 1 to 1^* will be considered a “jump.” Being $pij(tlag)=p(i→j|tlag)$ the conditional probability of a molecule in-state i to transition to state j after time t_lag, the lag-dependent transition matrix is defined as

Figure 9(b) presents a schematics of the resulting Markov chain. It is worth noticing that because of how 1^* is defined, a transition to or from 1^* by any other state except 1 is not possible, while $P1*,1(tlag)=1$⁠, regardless of the lag time.

The transition probability matrix is calculated for the case of $θ0=38°$ at equilibrium and for $Cacl=±0.02$⁠. The vector of time lags spans from 0 to 150 ps, sampled every 5 ps. The transition rates for jumping and rolling (that is $P1,1*$ and $P2,1$⁠, respectively) are shown in Fig. 10. We notice that contrary to what is usually observed in Markov models, the transition probability does not saturate, but rather starts to decrease after a maximum peak is reached. This behavior is explained by the use of a fixed control volume to track the molecules in the contact line region, which may exit the volume after a sufficiently long time and thus cannot be tracked anymore.

The transition probabilities are modeled with a simple exponential saturation curve $Pij=A[1−exp (−t/τ)]$⁠. A best least squares fit of the curve is performed for time lags below the peak of the transition probability, which is obtained by direct visual inspection and corresponds roughly to 60 ps for rolling transitions and 85 ps for jumping transitions. Albeit using large fit times lead to a slight overestimate of τ, they were taken this long to suppress noise. The estimates of τ and A allow us to both inspect the relevant timescale of molecular motion and the c.l. displacement rate. Results are shown in Table III. Roll events and jump events have two distinct time constants; the one of the former is comparable to the average lifetime of hydrogen bonds, whereas the one of the latter appears larger. The variability of time constants is not significant across capillary numbers, as all fall within the same range of uncertainty. The constants A show less relative variability and can be used to estimate P_ij for the value of t_lag corresponding to their maximum. Furthermore, we assume that the ratio of displacement rates would scale as the ratio of transition probabilities,

When estimating the exponential tilting using Eq. (8), we obtain a higher value for the rolling population than for the jumping one, as reported in the last row of Table III. Hence, we can speculate that indeed the displacement of the contact line occurs mostly due to rolling between liquid layers rather than jumping on the silica lattice. However, when comparing advancing and receding contact lines, the uncertainty is too large to determine if there is a relative difference in the rolling-to-jumping ratios,

Hence, we can state that while the two-layer structure is important in modulating the mobility of the contact line, the effect on its asymmetry is verified only indirectly by fitting the model to observables such as the dynamic contact angle and the contact line speed.

As stated in the Introduction, disagreement found in experimental investigations of spontaneous and forced wetting has promoted the idea that a local description of dynamic contact lines, i.e., a functional relation between U_cl and θ, cannot exist.²² Our simulation results clearly show that there exists an agreement across wetting modes and that contact line friction well models the local mobility of contact lines. In this section, we present some hypotheses on why experiments struggle to demystify dynamic wetting.

The optical resolution scale accessible by conventional microscopy techniques (which will be referred to as z_res) is usually reported to be in the range of 500 nm to 100 μm.^10,22–24 What can be directly measured in experiments is not the microscopic contact angle, but rather an apparent contact angle $θapp(zres)$⁠, which results from the interface bending at the resolution length scale. This contact angle, in general, differs from the microscopic one. According to the classical description of viscous interface bending, the mismatch may be quantified using Cox's formula,^5,6

with λ being a microscopic cutoff length, which is often identified in the slip length. Although attempts to disentangle the effect of viscous bending from the one of contact line friction using experimental observation have been tried,²⁵ the physical interpretation and the quantification of λ still remain difficult, especially for fluid-substrate combinations showing negligible slip. Hence, the disentanglement often has to rely on an arbitrary parameter.

It is reasonable to model the dynamics as overdamped and friction-dominated for droplets or menisci with linear dimensions of $∼10$ nm and for sufficiently small wetting speeds.²⁶ The same cannot always be said for the length scale accessible by experiments, where inertial-capillary effects may start to become relevant. Cox's theory is based on the assumption that the dynamics is viscous-dominated, i.e., small Reynolds number; the picture changes for moderate Reynolds numbers²⁷ and Eq. (10) loses its validity. Inertial under-damped dynamics is not self-similar; hence, experimental results are bound to differ unless the length scales, e.g., the droplet size, and the range of contact line velocities are both matched.

Because of the unavoidable limitations of experimental investigations, we wish to promote the idea of integrating molecular simulations in the validation process: an all-encompassing hydrodynamic model for moving contact lines is yet to be found and MD outputs could pave the way to its formulation. Already existing models can also greatly benefit from molecular benchmarks to perform parameter estimation, especially in the somewhat “pathological” case of no-slip surfaces. Following such intent, the results of the simulations conducted for this study are publicly accessible.^28–30 

Because of the resolution limit described above, it would seem impossible to obtain experimental evidence to distinguish between different contact line friction models. Nevertheless, if experiments are inadequate to directly characterize the relation between contact line speed and contact angle, it is still possible to quantify the contact line speed at which wetting failure occurs.³¹ Wetting failure may be defined as the condition for which a steady moving contact line cannot exist; it is realized either by the deposition of a thin film at the receding contact line or by vapor entrainment at the advancing contact line.³² 

The onset of wetting failure is typically identified by the condition $θapp=0$⁠, that is: the interface presents a zero-derivative inflection point at some distance from the contact line. It is possible to speculate whether the onset could be provoked by contact line friction. According to MKT, the maximum contact line speed is attained either for $θ=0°$ or $θ=180°$⁠; assuming that Eq. (10) is valid, it would imply $θapp=0°$ or $180°$ at any distance from the contact line. The two-layer model, on the other hand, predicts a maximum wetting speed attainable for some $θ*∈(0°,180°)$⁠. The condition can be obtained by simply setting

Figure 11 illustrates the contact line speed bounds consistent with Eq. (2). Assuming $θ*$ is far enough from 0° if the contact line is receding or 180° if the contact line is advancing, it would be, in principle, possible to observe no zero-derivative inflection of the interface at the onset of wetting failure. Under such conditions, even macroscopic measurements of θ_app could suffice to discriminate between contact line friction models.

In this work, we have expanded the application of the two-layer model by Johansson and Hess to cases of both advancing and receding contact lines, under spontaneous and forced wetting conditions, rationalizing the contact line asymmetry observed in Lācis et al. The two-layer structure of liquid water close to the three-phase contact line has been illustrated. Our results illustrate the importance of the disruption of the order and orientation of liquid molecules caused by the interaction with the solid wall: this effect needs to be accounted for in order to model contact line friction in situations such as the one of water wetting silica-like surfaces. Models that over-simplify local molecular conformation may be insufficient to describe wetting/de-wetting asymmetry. Baseline rates and probabilities for the rolling and hopping motions have been extracted from the transition time lags between the two identified conformations of water molecules in the contact line region. The results show the importance of the rolling of molecules in the second layers. However, unlike the previously studied case of rapid spreading,¹⁷ the sub-critical capillary numbers attainable on a hydrophilic surface are too small to obtain a good signal-to-noise ratio for the transition rates and distinguish the exponential tilting on advancing contact lines from the one of receding contact lines.

We have provided evidence to support the idea that the nanoscopic contact angle is a dynamic quantity and that its relation to contact line speed is invariant under the wetting mode in a regime where the dynamics is overdamped. The resolution accessible by atomistic molecular dynamics simulation allows us to quantify sub-continuous solid–liquid friction effects. It is important for experimental works incapable of resolving or inferring the curvature of liquid interfaces in the range of $∼10$ nm from the moving contact lines to acknowledge this limitation. Nevertheless, studying the onset of wetting failure could help distinguish between alternative contact line friction models even on the basis of macroscopic experimental observations. From the point of view of continuous fluid dynamics modeling of moving contact lines, allowing the microscopic contact angle to deviate from its equilibrium value is paramount to reproduce the correct contact line physics, unless $μ¯f≪1$⁠. Slip models that neglect contact line friction but manage to obtain quantitative agreement with MD or experiments should be regarded as effective, yet incomplete.

We finally comment on potential ways to expand on the current investigation. Even though the silica surface reproduces the fundamental inter-molecular interactions with water, surfaces in reality are not atomistically flat and homogeneous.³³ Topological and chemical defects or patterns can greatly influence contact line mobility; ways to adapt the contact line friction model to surfaces with defects or corrugations has been proposed in the past.^34,35 In this work, the contact line friction has been modified by tuning the interaction strength between the liquid and the substrate. Another interesting direction to explore would be to study the effect of changing the viscosity of the fluid.²³ Finally, it would be interesting to follow the thread of Sec. VI B and study how contact line friction affects the onset of wetting failure using molecular dynamics. This condition can be easily achieved by pushing the speed of solid walls beyond the critical droplet breaking value and observing contact line speed and interface curvature at the onset of instability.

See the supplementary material for details of the force field used in molecular simulation and for a thorough illustration of uncertainty quantification for contact line friction parameters.

We acknowledge funding from Swedish Research Council (INTERFACE Centre No. VR-2016–06119 Grant No. VR-2014–5680). We thank Dr. Petter Johansson for developing the code to perform data collection. We also thank Professor Shervin Bagheri, Professor Gustav Amberg, Professor Stephane Zalenski, Dr. Ugis Lacis, and Dr. Johan Sundin for the engaging scientific discussions which inspired the writing of this manuscript. Numerical simulations were performed on resources provided by the Swedish National Infrastructure for Computing (No. SNIC 2021/1–38) at PDC, Stockholm. Resources for data storage were provided by the Swedish National Infrastructure for Computing (No. SNIC 2022/23–46) and by Zenodo.

The authors have no conflicts to disclose.

Michele Pellegrino: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Berk Hess: Conceptualization (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The simulation output that supports the findings of this study is openly available in Zenodo at https://doi.org/10.5281/, Refs. 28–30. Analysis scripts in Python are available from the corresponding author upon reasonable request.