Nanoscale sheared droplet: Volume-of-Fluid, phase-field and no-slip molecular dynamics

The motion of the three-phase contact line between two immiscible fluids and a solid surface arises in a variety of wetting phenomena and technological applications. One challenge in continuum theory is the effective representation of molecular phenomena close to the contact line. Here, we characterize the molecular processes of the moving contact line to assess the accuracy of two different continuum two-phase models. Specifically, molecular dynamics (MD) simulations of a two-dimensional droplet between two moving plates are used to create reference data for different capillary numbers and contact angles. We use a simple-point-charge/extended (SPC/E) water model with particle-mesh Ewald electrostatics treatment. This model provides a very small slip and a more realistic representation of the molecular physics than Lennards-Jones models. The Cahn-Hilliard phase-field model and the Volume-of-Fluid model are calibrated against the drop displacement from MD reference data. It is demonstrated that the calibrated continuum models can accurately capture droplet displacement and droplet breakup for different capillary numbers and contact angles. However, we also observe differences between continuum and atomistic simulations in describing the transient and unsteady droplet behavior, in particular, close to dynamical wetting transitions. The molecular dynamics of the sheared droplet provide insight of the line friction experienced by the advancing and receding contact lines and evidence of large-scale temporal"stick-slip"like oscillations. The presented results will serve as a stepping stone towards developing accurate continuum models for nanoscale hydrodynamics.


Introduction
The motion of a two-phase interface over a solid surface has turned out to be a challenging problem for continuum fluid mechanics (CFM) models. This is most evident for models that assume a sharp interface between the phases and the no-slip velocity condition at the solid surface. Under these classical assumptions, one naturally ends up with an immobile line at which both fluid phases meet with the solid, so-called contact line. Clearly, a fixed contact line is in contradiction with, for example, our observations of a spreading drop on a surface, or of a liquid imbibition into a porous medium. This theoretical issue was identified already half a century ago by Huh & Scriven (1971) as a stress singularity at the contact line. Motivated by the importance of the moving contact line in applications such as printing (Kumar 2015), CO 2 storage and water management in fuel cells (Singh et al. 2019), a number of approaches have been suggested for overcoming the stress singularity in Navier-Stokes based solvers. The extensive work on the topic (Bonn et al. 2009;Snoeijer & Andreotti 2013;Sui et al. 2014) suggests that there exists no single continuum approach for describing the moving contact line. Instead, different models are suitable depending on the application and the representation of interfacial physics, and each model comes with its own sets of empirical parameters. Certain guidelines are required to determine these parameters.
Despite the vast amount of theoretical developments that have been presented (Voinov 1976;Cox 1986;Shikhmurzaev 1994Shikhmurzaev , 1997Kalliadasis et al. 2000;Eggers 2004;Flitton & King 2004;Wilson et al. 2006;Snoeijer 2006;Pismen & Eggers 2008;Pismen & Pomeau 2000;Snoeijer & Eggers 2010;Chan et al. 2012;Nold et al. 2018;Chan et al. 2020), a theoretical consensus is yet to be reached. To a large extent, progress is hindered by a lack of understanding the nanoscale physics of wetting phenomena (Afkhami et al. 2020;Afkhami 2021). For example, a fundamental question relates to the nanoscopic contact angle. The typical choice up till now has been a constant equilibrium angle (Kronbichler & Kreiss 2017). However, recent experimental (Deng et al. 2016) and numerical (Fernández-Toledano et al. 2021) evidence suggest a contact line velocity-dependent dynamic contact angle even at nanoscale. The development of more accurate measurement techniques is ongoing (Thormann et al. 2008;Eriksson et al. 2019), however, as of yet, we lack a complete insight into interface shape and velocity field near moving contact lines in nanoscale. Atomistic simulations have provided significant insight into the molecular physics at this scale (Gentner et al. 2003;Smith et al. 2016Smith et al. , 2018Perumanath et al. 2019;Lācis et al. 2020), although most systems used so far have been based on idealized force models between the liquid and the substrate.
In practice, there are few common approaches to numerically model moving contact lines in standard continuum methods. For sharp interface models, such as Volume-of-Fluid (VOF) and level-set (LS), the movement of the contact line is typically allowed by an explicit Navier-slip condition (Navier 1823;Spelt 2005) or by an implicit numerical slip at the contact line (Renardy et al. 2001;Afkhami et al. 2009). For diffuse interface models, such as the Cahn-Hilliard phase-field (PF) model (Jacqmin 2000), the contact line advances through diffusion even in the no-slip scenario. The Lattice-Boltzmann method (LBM) can be leveraged to solve different types of flow problems. The standard application of the LBM typically uses Shan-Chen (Shan & Chen 1993) or Gunstensen et al. (Gunstensen et al. 1991) models. It is argued that in the Shan-Chen model the contact line moves through phase change (Kamali et al. 2011), similar as in LBM with thermodynamically consistent potentials . The LBM method has been also used to solve phase-field equations , consequently, the motion of the contact line then occurs through diffusion. In the current work, we focus particularly on two of these models, namely, geometric VOF and Cahn-Hilliard PF models.
For the VOF method, different components to model the contact line (static, dynamic angle, hysteresis window, Navier slip, etc.) have been proposed over time. A recent comparison between different options can be found in Legendre & Maglio (2015). They concluded that the models incorporating dynamic contact angle better represent the experiments of a spreading drop. For the PF model, there are guidelines to select the model parameters through calibration with experiments (Yue & Feng 2011). A part of the calibration is choosing the diffuse interface thickness in a manner to satisfy the socalled sharp interface limit Magaletti et al. 2013;Xu et al. 2018). With this approach, a good agreement between PF simulations and capillary rise experiments has been demonstrated.
In general, however, the connection between the selected CFM models and molecular reality is not clear. It is not known how to choose the model parameters for an accurate representation of nanoscopic physics that determines the moving contact line speed. To address this, comparisons between molecular dynamics (MD) simulations and the chosen CFM models have been carried out (Qian et al. 2003; Barclay & Lukes 2019;Mohand et al. 2019). The typical MD work considers Lennard-Jones types of surfaces with large slip lengths. Only more recently, water MD simulations have been carried out by some of us (Johansson & Hess 2018;Lācis et al. 2020) over surfaces with negligible slip.
In this work, we generate benchmark data from MD simulations of wetting over noslip surfaces. These conditions can be reproduced employing SPC/E water on a smooth, silica-like substrate (Johansson et al. 2015). We choose a forced wetting set-up (Blake et al. 2015) -instead of a capillary-driven one (Villanueva & Amberg 2006) -due to the more versatile control of the wetting process. In particular, we choose a sheared droplet configuration (see figure 1a), which is a well-studied canonical problem (Jacqmin 2004;Sbragaglia et al. 2008;Gao & Lu 2013;Wang & McCarthy 2013) allowing simultaneous access to both receding and advancing contact lines. Furthermore, depending on the wall velocity U w , the system either i) stabilizes at a steady state (if U w < U w,c ) or ii) exhibits a non-trivial unsteady behaviour (if U w > U w,c ). Here, U w,c is the critical wall velocity describing the boundary between i) and ii). Hence, the sheared drop configuration provides rich interfacial dynamics that are challenging to capture with CFM models.
By assuming a negligible slippage, the space of input parameters is essentially reduced to wall velocities (U w ) and equilibrium contact angles between water and silica (θ 0 ). We then adopt a two-step approach. In the first step, we calibrate the continuum simulations against MD for a given pair (U w , θ 0 ). This yields the necessary PF and VOF parameters that best reproduce the steady droplet displacement measured from MD. In the second step, we fix the calibrated PF and VOF parameters and assess the predictive capability of the CFM models for different U w by characterizing the interface shape and the drop displacement both in the stable and unstable regimes. This is a extension of the work presented by Lācis et al. (2020), who evaluated CFM model performance in matching the MD results for a single steady wall velocity (U w < U w,c ) and a single equilibrium contact angle.
The molecular dynamics of the sheared droplet reported herein provide insight of the friction experienced by the advancing and receding contact lines. We demonstrate asymmetric features of advancing and receding lines and report evidence of large-scale temporal "stick-slip" like oscillations. These observations do not only enhance our physical understanding of moving contact lines, but also aid the development needed to increase the accuracy of continuum models.
The paper is organized as follows. In §2 we describe the flow configuration, CFM models and demonstrate the effect of the unknown parameters. The reference MD simulations of the sheared droplet system are described in §3. We calibrate the CFM models against MD in §4. Predictions from CFM are evaluated against MD in §5. Then, in §6 we provide insights into the molecular physics of the sheared droplet system. Following that, in §7 limitations of CFM models, fluid slippage, friction and potential future modelling directions are discussed. We conclude the paper in §8. In appendices A-E, important physical and technical details are provided.

Flow configuration and continuum models
We consider a two-dimensional system that is periodic in the streamwise direction and bounded in the vertical direction by two parallel horisontal plates located at y = 0 and y = H. A liquid drop of density ρ and viscosity µ is sandwiched between the plates such that its maximum width is W . The drop is surrounded by water vapour with density ρ v and viscosity µ v . The surface tension between the phases is constant and denoted by σ.
The numerical values of the geometry and fluid properties are reported in figure 1(a).
We study the response of the droplet to two parameters; i) U w , the constant velocity of the upper and lower walls moving in opposite directions and; ii) the equilibrium contact angle θ 0 . The former is interchangeably discussed in its non-dimensional form, using the Capillary number, The Capillary number, corresponding to critical wall velocity U w,c , is Ca c . For small Ca, the liquid slips on the solid, resulting in a steady deformed droplet shape as shown in figure 2(a). Above a certain critical Capillary number, Ca c , the liquid is deformed to such a degree that its interface to the surrounding vapour breaks, leaving behind multiple disconnected liquid droplets. The second control parameter, θ 0 , is a measure of the surface's affinity to water. In hydrophilic conditions (θ 0 < 90 • ), the affinity is strong, resulting in a larger droplet deformation compared to hydrophobic conditions (θ 0 > 90 • ). Under dynamic conditions, the contact line can be different from the equilibrium one. For advancing contact lines (marked with A in figure 2a), the liquid displace the vapour, while for receding contact lines (R in figure 2a), the vapour displace the liquid. Both these processes are largely determined by molecular interactions between the water and the substrate as depicted in figure 1(b).
The drop deformation induced by the moving walls is characterized with measures defined in figure 2(a). As a global measure, we introduce drop displacement of left (∆x l ) and right (∆x r ) two-phase interface, respectively. For more detailed characterization, we also evaluate the interface angle θ(y) for steady configurations. It is obtained from the slope of each linear segment on the interface.
In the continuum setting, the flow and pressure fields (u, p) are obtained by solving the incompressible Navier-Stokes equations in the domain containing both phases. In two dimensions, the equations are Here, f σ is the surface tension force, ρ (x, y) and µ (x, y) are spatially dependent density and viscosity, respectively. The functions ρ (x, y) and µ (x, y) take liquid and vapour values in the region occupied by each phase and undergo a sharp transition at the boundary between the phases. In this transition region, the volume force f σ is applied to model the force induced by the surface tension. Zero wall-normal velocity, u y = 0, is imposed on the moving impermeable walls. For the tangential velocity component, we impose a Navier-slip boundary condition where s is the slip length, u w is the wall velocity (U w at the top wall and −U w at the bottom wall) and n is the wall-normal coordinate. Although it is expected that s = 0 for the chosen liquid-solid combination, the implementation of the CFM models allows for tests with non-zero value. Equations (2.2-2.3) with corresponding boundary conditions are shared between the Cahn-Hilliard PF model and the geometric VOF model. However, the way surface tension and evolution of two-phase interface is treated differs.

Geometric Volume-of-Fluid model
In this model, a phase variable, C(x, y), is defined as 0 in the gas and 1 in the liquid. The phase variable thus represents the liquid volume fraction in each cell. It satisfies the convection equation ∂C ∂t = −u · ∇C. (2.5) The surface tension force is applied by the continuous surface force (CSF) method (Brackbill et al. 1992), where κ is the curvature of the interface. As a boundary condition for the phase variable, we impose a dynamic contact angle θ num . The contact angle is the angle between the interface (n i ) and wall (n) normals, and the sign might vary depending on convention. The relationship for θ num is inspired by Cox's theory, as described and evaluated by Legendre & Maglio (2015). It takes the form where G is a monotonically increasing function, λ is a microscopic cut-off length scale and ∆ is the wall-normal cell height. We confirm the grid convergence reported by Legendre & Maglio (2015) in appendix A.1. The capillary number Ca cl is based on the velocity of the contact line with respect to the wall. It is estimated by linearly interpolating the velocity   (Aniszewski et al. 2021). The full set of equations (2.2-2.3,2.5-2.6) together with the boundary conditions (2.4,2.7 and periodic inlet/outlet) are discretized on a regular cuboid grid with a staggered spatial representation. The cell spacing is constant in all directions. More details about the numerical method are given in appendix A.1. In order to compute the curvature near the wall for very hydrophobic droplets (θ 0 = 127 • ), we use a high-order scheme to approximate derivatives (appendix A.2 provides further details).
The only unknown free parameter for VOF simulations is the length scale λ. To test the effect of λ, we select θ 0 = 95 • and Ca = 0.20 < Ca c . The evolution of the drop displacement ∆x in time for λ values 0.47 nm, 0.66 nm, 0.94 nm, 1.87 nm, and 3.74 nm is shown in figure 2(b). Due to symmetry, ∆x l = ∆x r = ∆x in the CFM simulations, which we have verified numerically. We observe that larger λ correspond to smaller steady ∆x. As λ is increased, the contact angle at the surface (2.8) and thus the interface shape (2.7) near the wall is modified. The interface curvature is modified in such a way that the surface tension force across the interface opposes the friction force from the wall, which leads to a smaller ∆x. The ∆x for the initial time agrees between all λ values. In this short period, the contact line does not slip with respect to the substrate and the slope of ∆x(t) is equivalent to the wall separation velocity (2U w ).

Cahn-Hilliard phase-field model
We choose a model of a binary mixture with a classical fourth-order polynomial potential Ψ (see eq. B 1 in appendix B) (Jacqmin 2000;Carlson 2012). To describe the evolution of both phases, the PF model uses a phase variable C(x, y) ranging from 1 in the liquid to −1 in the vapour. At the interface, the function exhibits a smooth transition. The variable C is governed by a convection-diffusion equation (2.9) Here, M is the phase-field mobility and φ is the chemical potential. The latter is defined as (2.10) The chemical potential (2.10) contains the surface tension (σ) and the interface thickness ( ). The derivation of (2.9-2.10) is standard and can be found in Carlson (2012) and Jacqmin (2000). The surface tension σ is a physical parameter and is set to a value corresponding to the water liquid-vapour interface. The constants M and , on other hand, are typically treated as numerical parameters. Having solved for function C, the surface tension force in (2.2) is inserted as (2.11) The Cahn-Hilliard equation (2.9) is a fourth-order partial differential equation and thus two boundary conditions are needed on solid walls. The lowest order boundary condition is where µ f is contact-line friction and g is switch function (B 2). Equation (2.12) is also known as the non-equilibrium wetting condition and requires the equilibrium contact angle θ 0 as an input. Setting a non-zero µ f yields a dynamic contact angle different from θ 0 (Jacqmin 2000;Qian et al. 2003;Carlson et al. 2011). The second boundary condition on solid walls is ∇φ ·n = 0, which states that there is no diffusive flux through the walls.
The full system of fluid and phase-field equations (2.2-2.3,2.9-2.10) together with the boundary conditions is discretized and solved using open-source code FreeFEM (Hecht 2012). Adaptive mesh is used near the two-phase interface to capture the variation of the function C. More details about the PF equations and simulations can be found in appendix B.
The PF model has three unknown free parameters; the phase-field mobility M , the contact-line friction µ f , and the interface thickness . To provide an intuition about each parameter, we vary them one at a time. As before, θ 0 = 95 • and Ca = 0.20. Figure 3(a) shows ∆x for µ f = 0, = 0.7 nm and M = (3.5; 7.0; 10.5; 14.0; 17.5) × 10 −16 m 4 /(N s). We observe that the final steady ∆x value is reduced as M is increased, i.e. smaller deformation with increased diffusion. In figure 3(b), the displacement for = 0.7 nm, M = 1.75 × 10 −15 m 4 /(N s) and µ f = (0.0; 0.5; 1.0; 1.5; 2.0) µ is shown. Here, the ∆x increases with µ f and the contact-line friction has therefore an opposite effect compared to mobility. This competition has been previously clearly showcased by Yue & Feng (2011). Finally, the evolution of ∆x in time for = (0.18; 0.35; 0.70; 1.40; 2.80) nm is shown in figure 3(c) with µ f = 0 and M = 1.08 × 10 −15 m 4 /(N s). It can be observed that -in contrast to M and µ f variations -notable differences are present for initial times. As is decreased (in the direction against the arrow), the steady ∆x value is converging. This is a signature of the sharp interface limit Xu et al. 2018). For enriched understanding, in appendix F we report streamlines from PF near receding contact line for similar parameter variations.

Molecular dynamics simulations of the sheared droplet
We describe the polar molecules of the water droplet using the SPC/E model. This is the simplest model allowing hydrogen bonds with the solid substrate. It also offers an accurate description of water bulk and interfacial properties and retains a relatively low computational cost. The bounding walls are formed as mono-layers of SiO 2 quadrupoles (figure 4a) that are restrained into a hexagonal lattice. A quasi 2D system with depth 4.68 nm (figure 4c) is constructed. Albeit the composition of walls is structurally unrealistic, this simple surrogate configuration allows emulating the two fundamental electrostatic interactions characteristic to hydrophilic substrates (Johansson & Hess 2018). The first is the hydrogen bond between water and silica (dotted line in figure 4b). The second interaction is the adsorption of water molecules on the substrate. The adsorption occurs due to the attraction between water oxygen and silicon atoms (dashed line in figure 4b).
The strong electrostatic interaction is responsible for a very small hydrodynamic slip at the wall. Note that the current MD configuration does not include any chemical reactions that would occur at a real crystalline or amorphous silica surface.
The strength of the water-substrate interaction can be tuned by adjusting the charge distribution in SiO 2 (figure 4a). Different interaction strengths will result in different equilibrium contact angles θ 0 . We simulate the system via atomistic molecular dynamics in the NVT ensemble, using well-established force fields and thermostats. Details regarding the physical and numerical simulation setup can be found in appendix C.1.

Equilibration runs
First, we measure the equilibrium properties of the water drop between two static plates and generate thermodynamically consistent initial conditions. This is done through socalled equilibration runs (see appendix C.2). During the run-time of the MD simulation, we collect flow data in regular 0.2 nm × 0.2 nm bins for a time interval of 12.5 ps. This yields instantaneous density ρ i (x, y, t) and flow velocity u i x (x, y, t), u i y (x, y, t) data as functions of space and time. After the initial transient, ρ i (x, y, t) is averaged over time to reduce noise in the liquid-vapour interface shape. The interface shape is extracted based on this averaged density, i.e. ρ (x, y). The interface position is determined by seeking the location where the liquid density transitions from bulk density to very small vapour density. Example of density distributions are shown in figure 5(e,f). More details on interface extraction are provided in appendix C.3.
To determine the equilibrium contact angle, we extrapolate the interface angle towards a hydrodynamic wall position. We assume that the wall position is at the centre of the bin coloured in red (figure 5a-d). The q values are tuned to yield θ 0 = 127 • , 95 • , 69 • and 38 • . This allows us to investigate hydrophobic, neutral and hydrophilic wetting conditions. In parallel to θ 0 extraction, we also identify the first reliable bin for interface shape measurement. This bin is shown with green in figure 5(a-d). The interface angle computed from points closer to the wall exhibit extreme deviation from continuum description (figure C.2d). Consequently, for comparisons with CFM predictions (presented later), we always extract the interface shape from MD simulations neglecting the unreliable data points.
In previous investigations, density variations have been observed for Lennard-Jones (LJ) liquids near solid walls and two-phase interfaces (Bugel et al. 2011;Stephan et al. 2018). We determine the extent of ρ oscillations near the wall in our MD system. The water liquid density distribution along the height of the channel is computed as (3.1) Here, the boundaries for integration x l and x r are selected for each θ 0 to fall within the liquid phase for all y coordinates. The close-up near the bottom wall of ρ y (y) is shown in figure 5(a-d). Oscillations in ρ have smaller amplitude and occur over smaller distances than typically observed in LJ systems. The small layering is an outcome of the combination of the SPC/E water model and SiO 2 surface model. For the same value of θ 0 , the layering would propagate further into the liquid phase if an LJ surface would be used instead of SiO 2 .

Dynamic configuration
Having obtained the initial state from the equilibrium simulations, we turn to the dynamic configuration. Since the configuration is symmetric in a continuum sense, the final ∆x is obtained as an average between the left (∆x l ) and right (∆x r ) interface. We determine the drop displacement using the first reliable bins near the top and bottom walls (figure 5, green bins). Since there is no interface data near the hydrodynamic wall, the final steady drop displacement is obtained by extrapolating the interface shape. We use polynomial extrapolation as detailed in appendix C.4. Figure 6(a) shows ∆x(t) for θ 0 = 95 • for different capillary numbers, starting at Ca = 0.05 and then increased incrementally by ∆Ca = 0.05. For all Ca numbers up to and including Ca = 0.25, we observed a stable configuration. The obtained steady drop displacements for θ 0 = 95 • are summarized in second row of table 1. The table also reports displacement for the other equilibrium contact angles. Different Ca values are gradually tested until at least 3 simulations in a stable regime are gathered. For all θ 0 , as expected, we observe that ∆x increases with Ca.
The largest stable Ca = 0.25 at θ 0 = 95 • exhibits oscillations similar to the so-called "stick-slip" behaviour (Orejon et al. 2011;Varma et al. 2021). For a prolonged time, the contact line shows more resistance towards movement (i.e., stick). This period is followed by another in which the contact line exhibits less resistance towards movement (i.e., slip). We show a zoomed view of ∆x(t) for Ca = 0.25 in figure 6(c), where the stick-slip behaviour is identified with red arrows. Note that the partial stick-slip effect for (Ca, θ 0 ) = (0.25, 90 • ) is distinct from the oscillations observed for Ca = 0.20. We show the enlarged view of ∆x versus time for Ca = 0.20 in figure 6(b), where the red arrows depict the oscillation magnitude and time scale. We observe that the oscillations  at Ca = 0.20 are much smaller in magnitude and span smaller time scales compared to stick-slip-like motion (figure 6c). More discussion on the physics behind these oscillations is provided in §6.2. Similar oscillations with increased magnitude are observed also for θ 0 = 127 • at Ca = 0.90 and θ 0 = 38 • at Ca = 0.02. Simulations exhibiting stick-slip oscillations are denoted with * in table 1.

Droplet break-up
As the Ca number is increased further, ∆x(t) measured from MD does not stabilise around some finite value. Instead, ∆x(t) continuously grows. An example of ∆x(t) at (Ca, θ 0 ) = (0.30, 95 • ) is shown in figure 7. At 10 ns, the drop is only moderately deformed (figure 7a). At 23 ns (figure 7b), there are two drops at the top and bottom walls connected by thin thread of liquid water. Then, at around 23.5 ns break-up occurs and at 24 ns we observe two completely separated drops (figure 7c). Note that at the break-up instant, the slope of ∆x(t) changes distinctively. This is due to the absence of the surface tension force that was resisting the displacement. For the two separate drops, there is no competition between the friction at the top and bottom contact lines. Instead, the friction at the contact lines now ensures that the two drops follow the wall velocity and separate with speed 2U w . Hence, the critical Capillary number Ca c lies in between Ca = 0.25 and Ca = 0.30. In the sheared droplet configuration, the exact value of Ca c is determined by the most unstable contact line (either advancing or receding). For investigations of individual contact lines, other types of experiments should be performed, such as plunging/withdrawn plate Eggers (2005)

Calibrating CFM models against MD
The aim of this section is to identify the free parameters in PF and VOF such that the displacement ∆x obtained from the continuum models match the displacement obtained from MD simulations. The four calibration pairs (Ca, θ 0 ) are (marked with † in table 1) (0.60, 127 • ), (0.20, 95 • ), (0.10, 69 • ) and (0.015, 38 • ). For each θ 0 , we have chosen the largest steady Ca number available from MD simulations. If the chosen steady Ca number leads to large stick-slip like oscillations of ∆x (such as observed in figure 6c), we select the previous (smaller) steady Ca number. The unsteady simulations are not suitable for calibration, because ∆x(t) grows in time (figure 7).
Before calibration, we have to make sure that the same system -in terms of geometry and fluid properties -is represented in CFM and MD. The dimensions and physical properties are reported in figure 1(a). The bulk liquid density is obtained by taking an average of ρ y over the height of the channel This density is valid for all equilibrium angles. Viscosity and surface tension of liquid SPC/E water are taken from previous work (Lācis et al. 2020). The viscosity measurement of the vapour SPC/E phase from MD is impractical. Therefore, both viscosity and density of the vapour are determined from engineering tables (appendix D). For VOF simulations, we increase the vapour density to ρ V OF v = 9.9 kg/m 3 . This ensures numerical stability of the simulations, while keeping the influence of the vapour inertia negligible.
The size of the channel in the x-direction (distance between periodic boundary conditions) is matched with the distance between left and right periodic boundary conditions in MD system. The channel height is based on the chosen hydrodynamic wall position at the centre of the red bin (figure 5a-d), which results in H = 29.22 nm. This position  Table 2: Results of CFM calibration against MD. For each (Ca, θ0) pair, we report MD reference data (steady displacement ∆x) in third column. In following four columns, we report VOF parameters (eq. 2.4 and 2.8) and resulting steady displacements. In final five columns, PF parameters (eq. 2.4, 2.9 and 2.12) and resulting displacements are given.
is the subject of an investigation in itself (Herrero et al. 2019) and further discussed in appendix E. The slip length (2.4) used in CFM models in principle can be directly related to MD simulations. Previous work (Huang et al. 2008) has demonstrated that the slip length of MD system follows a quasi-universal relationship with respect to equilibrium contact angle θ 0 . However, for the chosen MD system (SPC/E water on SiO 2 surrogate wall), accurate slip length quantification has not yet been done. To obtain an indication about validity of s = 0 for the selected MD system, we compare streamwise velocity near the wall between PF and MD (appendix E). Through this comparison, we found that s = 0 holds for θ 0 = 38 • − 95 • , while for θ 0 = 127 • the appropriate choice is s = 0.44 nm. These values are used both in VOF and PF as input parameters without any additional fitting, see sixth and ninth columns of table 2. Note that s in VOF is fixed and distinct from λ -the length scale in condition (2.8), which is used for calibration.

VOF calibration
Figure 8 compares ∆x obtained from VOF (red) and MD (black) for the four calibration configurations. We have adjusted the parameter λ in (2.8) to find the best fit of ∆x as t → ∞. We observe that the droplet displacement from VOF stabilizes at values that agree well the atomistic simulations with the most challenging calibration configuration being θ 0 = 127 • (figure 8a). The λ values that reproduce the displacement obtained from MD are listed in sixth column of table 2. For the hydrophobic configurations, we obtained λ = 4.325 nm for θ 0 = 127 • and λ = 0.935 nm for θ 0 = 95 • . We observe that for θ 0 = 95 • , λ is roughly 4 times smaller than for θ 0 = 127 • . Naively, λ can be regarded as a slip-related length scale near the contact line. Larger λ for θ 0 = 127 • then suggests smaller friction, in line with the findings from PF calibration ( §4.2).
For θ 0 = 69 • and 38 • we found even smaller λ values, 0.313 nm and 0.146 nm, respectively. This continuous behaviour of decaying λ for smaller θ 0 , qualitatively, follows the contact line friction argument presented later ( §4.2). It is also interesting to note that the obtained λ values from this calibration procedure are of order 0.1 nm -4 nm. This is similar to what is used by Legendre & Maglio (2015), where they set λ = 1 nm in macroscopic simulations. The mesh spacing in the VOF simulations is ∆/2 = 0.457 nm. For θ 0 = 38 • and 69 • , we have λ < ∆/2 (table 2). On the other hand, for θ 0 = 95 • and 127 • , we have λ > ∆/2 (table 2). Here the sign of the logarithm in (2.8) is negative. This is typically not the case when applying Cox inspired relationships such as (2.8). Therefore imposed θ num for θ 0 = 95 • and 127 • should be viewed as numerical parameter to tweak the curvature near the wall.  Table 3: Calibrating PF with MD by changing µ f . Strategy proposed by Yue & Feng (2011).
In third column, we show the reference ∆x from MD. To the right, we show steady ∆x (in nm) obtained from PF for specified µ f /µ . If no steady state exists, we write "unst.". In all simulations, M = 3.5 × 10 −16 m 4 /(N s) and = 0.7 nm.

PF calibration
The droplet displacement ∆x produced by PF (green) and MD (black) are compared in figure 8 for the calibration configurations. The interface thickness, mobility and line friction that provide the best match with the displacement obtained from MD are listed in table 2. For PF, different combinations of parameters may provide a good fit; therefore it is prudent to have physically motivated guidelines.

PF calibration proposed in the literature
For the Cahn-Hilliard PF model, there is a standard calibration procedure proposed by Yue et al. (2010); Yue & Feng (2011). The sequential steps are; i) choosing the interface thickness that is suitable to describe the physical problem; ii) setting the PF mobility M according to the sharp interface limit  and; iii) calibrating the contact line friction µ f against experiments.
The interface thickness has to be smaller than the important physical length scales in the chosen system, which for the sheared droplet configuration are the water drop height H = 29.22 nm and width W ≈ 38 nm (figure 1a). Based on this, the interface thickness is set to = 0.7 nm. According to the sharp interface limit , the criterion for choosing M is Inserting the chosen values for interface thickness, vapour viscosity and liquid viscosity, we obtain M > 3.21 × 10 −16 m 4 /(N s). In the MD simulations, we do not observe any physical effect that would hint towards a large M value. Therefore, for all θ 0 in this section, we select M = 3.5 × 10 −16 m 4 /(N s), which is close to the lower limit. The last step is to carry out simulations with different µ f until a steady state is attained that matches the displacement obtained from MD. Table 3 summarizes the displacement obtained for each calibration couple (Ca, θ 0 ) at µ f = 0, 0.1, 0.2, 0.3, 2.36 and 11.84. The fourth column of table 3 shows the displacements for µ f = 0. We observe that no steady solution has been obtained for θ 0 = 127 • . The steady ∆x obtained for θ 0 = 95 • overestimates the MD value (third column of table 3) by roughly 50%. Increasing µ f only deteriorates the agreement. Therefore, we conclude that the matching procedure proposed by Yue & Feng (2011) is not adequate to calibrate the PF model for the chosen nanoscale configuration at θ 0 = 127 • and 95 • .
For θ 0 = 69 • and 38 • , on the other hand, µ f = 0 results in an underestimated drop displacement (table 3) compared to MD. We observe that the required PF contact line friction for θ 0 = 38 • (µ f = 11.84 µ ) is roughly five times larger than the one for θ 0 = 69 • (µ f = 2.36 µ ). More hydrophilic θ 0 entails a stronger affinity between the liquid and the wall. Stronger affinity, in turn, can yield larger friction, which is consistent with the obtained µ f values.

4.2.2.
Calibrating PF by adjusting mobility for θ 0 = 127 • and 95 • A way to reduce friction near the contact line is to increase M (figure 3a) and thus allow for more diffusion. Therefore, we impose µ f = 0 and increase M for θ 0 = 127 • and 95 • until the steady ∆x from the PF agrees with MD. The required PF mobility values are reported in table 2. We observe that for θ 0 = 95 • , a three times larger M is required compared to the M set by the sharp interface limit (see tenth column of table 2). Whereas for θ 0 = 127 • the M value has to be increased by a factor of hundred.
Note that this calibration procedure is compatible with the sharp interface limit. The condition (4.2) can be rewritten in terms of interface thickness as < 4 M 1/2 µ 1/4 v µ 1/4 . As we increased M , the interface thickness = 0.7 nm was kept constant, and therefore the condition is satisfied. This calibration procedure produces however a noticeable diffusion near the contact line (figure F.1c), which is not observed in the MD simulations.

Predictions from PF and VOF models
We have shown that continuum systems can be tuned to match the final steady droplet displacement computed from MD simulations. It is also necessary to understand how well the CFM models capture other key features of the system, including the interface shape and the time-dependent transient behavior of the droplet. Moreover, an important practical aspect is the accuracy of the CFM when it comes to predicting the droplet behavior away from calibration conditions. The aim of this section is therefore to quantitatively characterize the sheared droplet system for a range of capillary numbers. Specifically, in this section, we fix the parameters for PF and VOF to values reported in table 2. Figure 9 shows the droplet displacement as a function of time for a fixed θ 0 = 95 • but different Ca. Figure 9(d) is identical to figure 8(b) and corresponds to the conditions for which the system was calibrated for, i.e. Ca = 0.20. We observe that both PF and VOF capture the transient dynamics very well overall. The PF model is slightly slower (predicts smaller ∆x at the same time instant) than VOF and MD. Figure 9(a-c,e) shows the transient behavior of CFM models in off-calibrations conditions for Ca = 0.05, 0.10, 0.15 and 0.25. For moderate Capillary numbers Ca 0.10, we observe qualitatively similar behaviour as for Ca = 0.20. The CFM models predict the transient and steady ∆x values rather accurately. The PF, however, is always slightly slower compared to the VOF model. For Ca = 0.05, the steady ∆x value is slightly lower in CFM models than in MD. The PF model is slower in the transient compared to the VOF, and the agreement with MD is arguably worse.

Time evolution of drop displacement
To conclude θ 0 = 95 • investigations, we consider the CFM model predictions for the unsteady configuration with Ca = 0.30 (figure 9f). The results from PF, VOF, and MD are indistinguishable for the first 5 ns showing excellent predictive capability. For later times VOF simulation over-predicts and PF simulation under-predicts the ∆x observed from MD. This is in line with observations in steady situations (figure 9a-e). Both VOF and PF exhibit the rapid change of slope at around 23.5 ns, corresponding to the drop break-up (discussed in §3).
We have carried out the same investigation for θ 0 = 38 • , 95 • and 127 • . Qualitatively similar results to those observed in figure 9 are obtained, see Supplementary Figures 1-3.

Interface shape
For interface shape comparisons, the data is presented as the variation of the angle along the interface θ (y) (figure 2a). The MD data in Figure 10 does not extend all the way to the wall, so that the exact value of the dynamic contact angle is not known. At θ 0 = 95 • , the dynamic contact angle in PF simulations is equal to the equilibrium angle since µ f = 0 (table 2). For VOF, the dynamic contact angle is different than θ 0 . However, the angle is larger (smaller) than the equilibrium angle for the receding (advancing) contact line. This is opposite to the understanding arising from analysis of uncompensated Young stress. The source of this effect is the value of λ = 0.935 nm > 0.467 nm = ∆/2 (table 2).
We have repeated this comparison of interface shapes at different Ca numbers between MD, PF, and VOF for equilibrium contact angles θ 0 = 127 • , 69 • , and 38 • . The results are reported in Supplementary Figures 4-6 and are similar to what is presented above. Figure 11 shows the steady displacement as function of Ca for different θ 0 . As expected, the MD, PF and VOF points collapse for the calibration configurations (marked with vertical arrows). For θ 0 = 127 • configuration, the PF and VOF model overestimates the steady displacement for large Ca. A possible reason behind the discrepancy between PF, VOF and MD could be the slip length s . It has been determined for Ca = 0.60 and kept constant for runs with different Ca. In general, the slip length can increase for larger wall shear stress (Thompson & Troian 1997). This, consequently, would reduce ∆x and possibly move PF and VOF predictions closer to MD results. For θ 0 = 38 • , the agreement between PF, VOF, and MD diverge as we increase Ca > 0.02. A reason of the disagreement could be the increased contact line friction ( §6.1) at receding contact line, from which liquid film is formed ( §5.5 and figure 14). This asymmetry is not taken into account in the CFM simulations.  Figure 12 shows the critical Capillary number, Ca c , as a function of θ 0 from the three methods. To determine Ca c for a given θ 0 , we take the mean of the largest steady (Ca s ) and smallest unsteady (Ca u ) Capillary number that was simulated, i.e.

Steady drop displacement
The uncertainty is determined as half of the difference between these Ca numbers.
To reduce ∆Ca c in CFM, we sample the Ca space with smaller intervals. It can be observed from figure 12 that Ca c increases for larger θ 0 , which has also been reported in previous numerical (Sbragaglia et al. 2008) and analytical (Hocking 2001;Eggers 2004) investigations. We did not manage to run any VOF simulations for θ 0 = 127 • and Ca 0.80. Due to a large λ (table 2), the dynamic angle for Ca 0.80 falls outside of physically admissible range (0 • , 180 • ). Consequently, in figure 12 we compare only PF prediction with the MD data at θ 0 = 127 • . The PF model slightly under-predicts the Ca c value. This discrepancy could again be due to the uncertainty in the slip length s . Figure 12 shows that predictions of Ca c for 95 • and 69 • agree with MD results within the accuracy bounds of the MD. For θ 0 = 38 • , however, discrepancies are observed, where Ca c computed from both CFM models are larger compared to the MD results. Note that the differences at 38 • are enhanced due to the logarithmic scale of Ca c axis in figure 12.

Above the critical Capillary number
In this section, we investigate the accuracy of the CFM models for predicting the unsteady breakage of the droplet. This constitutes the most challenging test of the CFM models in off-calibration conditions. Two configurations are selected with Ca > Ca c , namely, θ 0 = 95 • at Ca = 0.30 and θ 0 = 38 • at Ca = 0.05. Figure 13 shows the drop shape at four time instances for (Ca, θ 0 ) = (0.30, 95 • ). We observe that the three models provide similar deformed states at t ≈ 5 ns and t ≈ 10 ns, see top two rows in figure 13. The time instance just before and right after the break-up is shown in third and fourth rows in figure 13. The thread connecting the lower and upper drop is very thin in MD and VOF simulations, compared to the PF simulation. In addition, the thread in VOF is comparably long and exhibits grid-to-grid like oscillations. There are also pronounced differences in the neck of each satellite drop -the region in which the thread transitions to the drop shape. The PF simulations show the thickest neck (figure 13g), followed by MD with slightly thinner neck (figure 13c) and VOF with

Molecular physics of the sheared droplet
In this section, we present molecular phenomena of the sheared droplet that are particularly challenging to model in a continuum model. First, we present contact line friction measurements directly from MD and asymmetry between advancing and receding lines for hydrophilic and hydrophobic configurations. Second, we discuss the nature of the stick-slip like oscillations. Finally, we show the inevitable three-dimensional nature of the drop break-up.

Contact line friction measurements from MD
To extract the µ f from MD, we use the formula proposed by Yue & Feng (2011), Here, θ is the dynamic contact angle at the wall. This expression is an approximation of the wetting boundary condition (2.12) in case of no-slip and small Capillary number. Note that there are alternative approaches to determine contact line friction, for example, based on equilibrium simulations, as proposed by Fernández-Toledano et al. (2019. For this work, however, we have determined that non-equilibrium approach based on fitting (6.1) to MD data is the most efficient approach.
To determine the µ f from (6.1), the dynamic contact angle θ has to first be determined  for each Ca number. We extract the dynamic contact angle above the first reliable bin (green bins in figure 5a-d). To reduce the noise, we use polynomial interpolation of the interface shape to read the dynamic contact angle at the chosen location. For consistency, we also re-evaluate the θ 0 at the same distance from the wall. The obtained dynamic contact angles for all θ 0 values are gathered in figure 15. For each θ, we also display error bars, obtained by adding ±2 to the polynomial order of the interpolation. As expected, when Ca number increases so does the deviation of θ from θ 0 . We use least-squares fit to match (6.1) to θ measurements. The fit provides µ f /µ values and error intervals. We fit the measurements taken at the top-left/bottom-right and top-right/bottom-left contact lines separately. This allows the observation of an asymmetric line friction. The obtained best-fit lines for all equilibrium contact lines are shown in figure 15. The contact-line friction are listed in table 4 along with previously reported µ f from calibration of PF simulations.
We observe that larger line friction parameters are measured for smaller θ 0 . This observation is consistent with the molecular-kinetic-theory (MKT). It states that the line friction scales with the work of adhesion needed to desorb water molecules from the substrate, which in turn increases as the surface becomes more hydrophilic (Blake & Haynes 1969). Interestingly, we observe a difference between advancing and receding line friction for θ 0 = 127 • and 38 • . Future and dedicated work is required to investigate this asymmetry in depth. In particular, to determine if the asymmetry depends on the position of the hydrodynamic wall. The observation of asymmetry is not specific to (6.1). A linearised MKT model -where line friction is defined in a slightly different waywould not change the conclusions.
The line friction obtained by fitting expression (6.1) directly to MD data are larger than the ones obtained through calibration of PF against MD (table 4). This fact seems to entail some missing physical effects in the PF model. Ideally -in a potentially more advanced PF model -one would employ µ f obtained from MD directly in the PF boundary condition (2.12).

Stick-slip-like oscillations
To differentiate between fluctuations caused by molecular-scale motion and large stickslip-like oscillations we define three reference length scales. The characteristic length scale of interface fluctuations far from contact lines can be estimated by balancing the thermal energy k B T with the energy due to surface tension σl 2 th , giving l th = k B T /σ 0.27 nm. The other two scales are the van der Walls diamater of SPC/E water molecules σ SP C/E 0.32 nm and the hexagonal lattice spacing of the substrate d hex = 0.45 nm. We then compute the root-mean-square for fluctuations of the contact line displacement RMSF= ( ∆x 2 − ∆x 2 ) 1/2 for each of the four contact lines, both at equilibrium and We observe that for Ca < 0.25 the RMSF is comparable with d hex , hinting that the observed fluctuations for moderate Capillary numbers are caused by the same process producing fluctuations at the equilibrium, that is the local thermal-induced pinning (depinning) of the contact line on (from) adsorption sites close to the average contact line position. It is worth noticing that the scale of these fluctuations is larger than the one expected on the interface far from contact lines. This observation is confirmed by examining the RMSF across the whole interface at equilibrium (see Supplementary Figures 7-10). Note that for some contact angles, the fluctuations in equilibrium simulations are smaller compared to d hex .
We examine closer the MD simulation with θ 0 = 95 • at Ca = 0.25 (figure 6a), which shows a much larger contact line RMSF, incompatible with lattice spacing driven fluctuations. The speed of the drop displacement during the stick-slip like motion is much smaller than U w (figure 16b). We conclude that the contact line does not entirely pin when resisting motion and only partially slips when complying with it. Moreover, for most of the time evolution, ∆x l and ∆x r are synchronized in an oscillatory motion. However, there are a few time intervals (between 17 ns and 20 ns; 45 ns and 47 ns) where indeed complete stick-slip occurs. In these intervals, the advancing and receding speeds match the magnitude of wall velocity, see inset of figure 16.
It appears that a local (pinning/depinning) and a global (oscillations) processes coexit. Pinning/de-pinning can be explained by the fact that Ca is close to the critical Capillary number, Ca c . The physical interpretation of the global oscillations is more delicate. In our modelling approach, we have implicitly assumed that the contact line motion is over-damped. This means that there should be a (possibly nonlinear) direct relation between force and speed, in which the proportionality constant is the contact line friction. This may not be true for large wall velocities. In such a case, the effects of the neighbouring flow inertia come into play. In this scenario, the coupling between positions and forces is more complex. The stochastic forcing produced by thermal fluctuations of the microscopic contact angle is no more completely dissipated. Instead, these may excite oscillation modes of the whole interface. Stick-slip has been theorized for flat surfaces and homogeneous fluids under some flow conditions (Hocking 2001;Eggers 2005;Varma et al. 2021). However, to the best of the author's knowledge, it has not been directly observed.
The selected standard CFM models are not capable of describing stick-slip like oscillations. The Navier-Stokes equations, underlying the PF and VOF models, are inherently deterministic, where all the thermal oscillations are averaged out. To model stick-slip like processes in a continuum model, a possible approach could be to introduce random fluctuations on the imposed θ 0 . The distribution of the contact angle oscillations has been previously identified (Smith et al. 2016). Fluctuations of the contact angle would in turn induce oscillations in ∆x(t).

Three-dimensional nature of drop break-up
In figure 13(d,h), we have observed that the interface shape obtained from MD and PF display more diffuse regions at the tip of the broken thread. It is tempting to conclude that PF correctly captures the MD behaviour. However, the MD snapshots in figure 13 have been averaged in spanwise z−direction, while PF simulation is a pure 2D result.
To obtain a better physical understanding of the exact break-up process, we investigate the molecular field of the MD system exactly at the break-up. In figure 17 we show water molecule locations for θ 0 = 95 • at Ca = 0.4 shortly before the break-up. Observing the drop from the side (figure 17a) it seems that the thread is not interrupted. In reality, the thin thread develops three-dimensional holes (figure 17b) before disconnecting completely. Averaging in the z direction results in a lower density (than the bulk) and thus giving the impression of a diffuse interface (figure 13e).
It is also worth noticing that the formation of these three-dimensional threads occurs very quickly. The time it takes for the threads to form is around 50 ps before the actual break-up. The break-up itself occurs after a few nanoseconds from the instant the twodimensional neck starts to form. We can thus infer that for the chosen combination of a molecular system and domain size there exists a time scale separation between 2D and 3D breakage dynamics. This in turn suggests that the selected MD system dimensions are reasonable to produce 2D results and can be directly compared with 2D CFM simulations.

Discussion
Based on the data presented in the previous sections, we discuss the accuracy, limitations and future development directions of CFM models when it comes to modeling the molecular physics in a nanoscale channel.

Challenges of the chosen CFM models
For θ 0 = 127 • , we observed differences in the steady drop displacement for Ca > 0.6, see figure 11(a). In general, very high contact angles are challenging to model using the particular CFM parameters we have chosen for calibration. For VOF, it is important to consider condition (2.8) as a numerical model that allows adjusting the interface curvature at the contact line. For θ 0 = 127 • this was achieved by increasing λ to values larger than the grid size, i.e. λ > ∆/2 (table 2). This leads to a negative sign of logarithm in condition (2.8). This means that for the bottom-left receding contact line (figure 2a), where Ca cl is negative, we have G(θ num ) > G(θ 0 ) and thus θ num > θ 0 . In other words, the imposed numerical contact angle is increased such that the curvature at the receding contact line results in a sufficiently large positive forcing in the x-direction on the fluid. The imposed numerical contact angle does therefore not approximate the true dynamical contact angle of the system. For PF, the mobility parameter M is calibrated ( §4.2.2) to increase contact line velocity for high θ 0 . Here, mobility is considered as numerical tuning parameter to match the droplet displacement and does not correspond to the actual molecular diffusion at the interface (see appendix F and figure F.1a-c).
Very low contact angles impose other challenges for CFM. For θ 0 = 38 • , we observed that steady drop displacement diverged (figure 11d) between the three models. Indeed, the prediction of critical capillary number showed deviations between the methods (figure 12). In addition, the time it takes for the drop to evolve in different shapes (figure 14) is different in all simulation methods. This is despite the fact that the λ obtained through VOF calibration ( §4.1) is reasonable (λ < ∆/2) and the standard PF calibration procedure ( §4.2.1) works well. By investigating MD data directly, asymmetric contact line friction was observed (table 4). For θ 0 = 38 • we observed significantly larger µ f for receding contact line compared to one for the advancing contact line. In the CFM models, on the other hand, the contact line properties (µ f for PF and λ for VOF) were the same for both advancing and receding sides. Since the receding contact line becomes unstable first (figure 14), this is the likely reason of the CFM model inaccuracy.
Finally, there are limitations of the chosen CFM at very low and high Ca numbers. At high Ca numbers (close to Ca c ) we observed enhanced oscillations of ∆x(t) (figures 6c and 16, stars in table 1). More detailed analysis of these oscillations were presented in §6.2 for θ 0 = 95 • . The CFM models does not include the intrinsic oscillations present in the molecular reality. Consequently, the agreement between the CFM and MD in drop time evolution (figure 9e and Supplementary Figure 6) is degraded. This could also be a potential source of increased discrepancy between interface shapes from CFM models and MD (figure 10e). Also for low Ca numbers (for example, Ca = 0.05 in figure 9a), we have observed relatively large difference between CFM models and MD results. The underlying cause for this inaccuracy is the relatively small drop displacement. As seen in figure 6(a), the MD oscillations have roughly the same magnitude for all Ca numbers.
Consequently, the signal to noise ratio in MD is much larger for smaller Ca numbers and the large oscillations can give impression of larger inaccuracy of CFM models.

Fluid slippage and contact line friction
In this work, the nanoscale molecular system has a negligible hydrodynamic slip. It was observed that the slip length s exhibits only small variations to θ 0 . Below bulk liquid, the slip length had to be 0.44 nm for θ 0 = 127 • and 0 nm for all other θ 0 . As discussed in §7.1, for low-friction configurations θ 0 > 90 • , it was necessary to adjust λ and M to further enhance contact line movement relative to the wall. This serves as a hint that the friction near the contact line is much smaller compared to what is modelled through bulk slip length s . In general however, other physical mechanisms that are intrinsic to the contact line are expected to co-exist with molecular slippage.
We have extracted contact line friction directly from MD (table 4). For hydrophilic substrates, a particularly high friction value was obtained. This is consistent with the formation of a microscopic water film for a relatively small capillary number (i.e. θ rec ∼ 0 for Ca 1). On the other hand, the interpretation for the almost vanishing receding friction on the θ 0 = 127 • surface is less obvious. To the best of the author's knowledge, nowhere before an asymmetric behaviour has been reported. It is tempting to explain the asymmetry by stating that hydrophilic surfaces are easier to wet rather than de-wet and vice-versa for hydrophobic surfaces. This puts the classical view -that hydrophilic substrates are high-friction surfaces and that hydrophobic substrates are low-friction surfaces -under doubt. Indeed, it is not clear whether the value of contact line friction can be predicted from θ 0 alone (Liu et al. 2019b;Wang 2019). Reasoning with the frame of mind of molecular-kinetic-theory instead, line friction asymmetry hints toward some complex physics modulating adsorption/desorption of molecules at the contact line. Fluid/surface interface energy alone is not sufficient to describe asymmetry between adsorption and desorption. More in-depth examination of sub-continuum fluid displacement near contact lines will be required to arrive at a physical understanding of this process. When the asymmetry is understood, it is straightforward to impose different µ f values in PF and different λ values in VOF for advancing and receding contact lines.

Potential modelling directions
It has been previously proposed that a better way to model the moving contact line is to use the so-called generalized Navier boundary condition (GNBC) as first hinted by Blake (1993). This approach has been later on evaluated against MD simulations (Qian et al. 2003(Qian et al. , 2006Mohand et al. 2019) and good agreement has been found. However, recently Lācis et al. (2020) have tried to match GNBC with MD simulations exhibiting negligible slip but were not successful in demonstrating any advantage of GNBC against no-slip and Navier-slip boundary conditions.
In this work, we also use substrate with negligible hydrodynamic slip ( s = 0.44 nm for θ 0 = 127 • and s = 0.0 nm for θ 0 = 38 • -95 • ). Previously (Lācis et al. 2020), for θ 0 = 95 • we extracted s = 0.17 nm, which is close to what we have in the current work. However, the most important overlap between the current and our previous (Lācis et al. 2020) work is that the slip length imposed at the solid wall is constant over the surface. This is the main reason why the GNBC for the selected system (Lācis et al. 2020) did not exhibit any advantage. Very small slip length corresponds to very large friction at the contact line. Consequently, the addition of uncompensated Young's stress does not lead to significantly modified flow near the contact line.
Nevertheless, there are no solid arguments as to why the effective friction exactly at the contact line should be the same as below the bulk liquid. It could very well be that the effective slippage (friction) exactly at the contact line must be prescribed larger (smaller) compared to below the bulk liquid. This has the potential for improving both PF and VOF ability to match the MD results. An alternative approach to prescribing larger friction in VOF simulations would be to use so-called staggered slip or negative slip (Hartmann et al. 2021). Another effect, which was not considered in this work, is so called disjoining pressure (Pismen & Pomeau 2000). Including it would also allow for direct modification of contact line motion. Furthermore, more detailed studies of local rheological effects (including orientation parameter of water) could provide insight of detailed mechanisms governing film formation and contact line friction. An accurate description of Navier slip related friction near the contact line could improve results attainable also using the GNBC condition. This is because the friction parameter is an important input in the GNBC. Infinite friction parameter renders GNBC ineffective, whereas gradually reducing friction parameter (increasing slip) would amplify the GNBC effect on the velocity near the contact line.
Another aspect is that the MD results hint towards asymmetric properties of advancing and receding contact lines. Consequently, contact line friction and possibly local slip length could be different. This enlarges the parameter space enormously. Fundamental investigations into the possible cause of asymmetry between the adsorption and desorption process would be welcome. These could potentially shed more light on what input should be given to CFM models to match the molecular reality. Alternatively some kind of hybrid methods that allow the matching between MD and Navier-Stokes solvers (Hadjiconstantinou 1999;Zhang et al. 2017;Borg et al. 2018) could be used. These approaches would alleviate the need to understand the asymmetric properties of the contact line and provide direct coupling between MD and CFM solvers.
Finally, there are other CFM models available that could be benchmarked against the molecular data produced through this work. For example, there exist different PF models, such as van der Walls (Laurila et al. 2012) or Cahn-Allen (Eggleston et al. 2001). The level-set (Tornberg & Engquist 2000) model or Lattice Boltzman (Chen et al. 2014) method are other potential candidates for the simulation of multiphase flows; Latva-Kokko and Rothman (Latva-Kokko & Rothman 2007) showed for instance how a no-slip LB model is able to automatically capture the speed-dependent dynamic contact angle and the interface-local slip length. The number of freely adjustable parameters differs between all alternatives. However, the issue of the appropriate velocity boundary condition will be shared between all of models based on single continuum description. A good recent classification of multiphase models can be found in work by Soligo et al. (2021). It has to be recognised that the hybrid models (Zhang et al. 2017;Liu et al. 2021) alleviate the need to understand the fundamental mechanisms near the moving contact line and provide a way to couple MD and CFM directly. This is another alternative that should be considered for efficient simulations of multiphase systems. If the stickslip like oscillations of the contact line are important, some other means of continuum modelling can be pursued. For example, the fluctuating hydrodynamic interfaces model proposed by Flekkøy & Rothman (1996);Smith et al. (2016) could be evaluated as a suitable choice. Lastly, we mention the possibility of regularizing the contact line stress singularity via the Brinkman's model for porous surfaces (Devauchelle et al. 2007).
Despite the missing physical mechanisms in the PF and the VOF, we have demonstrated that sufficiently accurate predictions of interface shape, drop displacement and critical Capillary number can be obtained. The only prerequisite is that the PF and VOF simulations have to be calibrated with the MD once for given θ 0 .

Conclusions
We have calibrated a standard Cahn-Hilliard phase-field model, as well as a standard geometric Volume-of-Fluid model with Cox-like wetting condition, against molecular dynamics simulations of water over a no-slip substrate. The no-slip behaviour in the MD system is an outcome of electro-static bonds between polar water molecules and polar wall molecules. Two parameters (mobility M and contact line friction µ f ) were adjusted in PF simulations and one parameter (Cox cut-off length scale λ) in VOF simulations. Four different equilibrium contact angles (θ 0 = 127 • , 95 • , 69 • and 38 • ) were investigated. For each θ 0 , the largest steady stick-slip free simulation was selected for calibration. The PF and VOF models were calibrated to match the steady ∆x -single scalar macroscopic measurement -observed in MD. Using the calibrated parameters, a series of simulations were carried out for other solid wall velocities. We demonstrated that CFM simulations can sufficiently accurately predict the drop displacement without any additional adjustments. The critical Capillary number predictions from CFM models also displayed good agreement with MD data.
In addition, we have showcased predictions for two unsteady sheared droplet configurations of θ 0 = 95 • and θ 0 = 38 • . The CFM models predicted all qualitative features of the MD simulations. For θ 0 = 95 • , drop displacement and break-up in half were predicted. For θ 0 = 38 • , thin film deposition and coalescence with the periodic image were predicted. Despite the quantitative differences in the time instances of the similar shapes, the CFM predictions exhibited good agreement with MD results.
We identified molecular physics that to the best of the authors' knowledge have not been previously reported. We extracted line friction directly from MD and compared it with the PF calibration results. The MD results showed the same trend as obtained with PF. In addition, the resulting contact line friction µ f values were asymmetric between advancing and receding contact lines for θ 0 = 127 • and 38 • .
Finally, we have discussed the variations of PF and VOF parameters for matching the MD results for all θ 0 values. We identified that the currently chosen CFM models seem to be lacking a way to describe reduced friction near the contact line for θ 0 95 • . A possible future direction to remedy this shortcoming would be to introduce larger slippage (lower friction) locally near the contact line. In a continuum setting, this could correspond to having a spatially varying slip. In addition, we anticipate that the asymmetric behaviour of the advancing and the receding contact lines is the source of the inaccuracies in PF and VOF when predicting the θ 0 = 38 • results.
We have demonstrated that by calibrating the CFM once for each θ 0 by targeting a single global measure ∆x it is possible to obtain many accurate predictions of interface shape, ∆x as a function of Ca and also prediction of Ca c . This is despite the fact, that the selected MD configuration exhibits practically negligible slippage. The selected CFM models seem to perform very well for close to neutral and slightly hydrophobic configurations (θ 0 = 95 • and 69 • ). On the other hand, more deviations were observed for hydrophobic (θ 0 = 127 • ) and hydrophilic (θ 0 = 38 • ) configurations. These accuracy limits have to be kept in mind, if these CFM models are applied in similar conditions.
The results of this study continue to enrich our understanding of connections between continuum mechanics simulations and molecular reality. We believe that this work provides important insights into PF and VOF models, associated open questions, and the required calibration procedures. Properly calibrated, both PF and VOF can serve as useful tools for investigations of technological applications.

Supplementary data
Supplementary figures, movies and data files for easier figure reproductions are available at ... analysed data. U.L., M.P., J.S. and S.B. wrote the paper with feedback from all coauthors.

Appendix A. Details of the geometric Volume-of-Fluid model
In this appendix, we provide additional details of the VOF model. First, we describe the numerical implementation of the solver. After that, customisation for the θ 0 = 127 • configuration is explained.

A.1. Numerical implementation
All VOF simulations used a resolution of N x = 256 (∆ x = 0.624 nm) cells in the streamwise and N y = 32 (∆ = ∆ y = 0.913 nm) in the wall-normal direction. PARIS solves the general three-dimensional equations. The two-dimensional behaviour was obtained using a thin domain, two cells wide in the spanwise direction (with periodic boundary conditions). We performed the simulations using a first-order time scheme, and the pressure was computed using the HYPRE library. Momentum was advected with a second-order central difference scheme. Equation (2.5) was solved using the builtin implementation of the algorithm by Weymouth & Yue (2010), computing the fluxes of C on the faces of each cell and updating C accordingly. The equation for the dynamic contact angle (2.8) was solved using the implementation reported by Sundin et al. (2021). Boundary conditions were implemented through a ghost layer (a layer of cells outside the computational domain where numerical quantities can be imposed).
The curvature of the interface was calculated using height functions. Height functions give the distance to the interface from a reference plane aligned with the grid. The values of the height function in specific cells, called heights, are computed by integrating C. Wall-parallel (wall-normal) heights provide interface x-coordinates (y-coordinates) for equidistant y-locations (x-locations) corresponding to the simulation mesh. The dynamic contact angle was imposed by prescribing wall-parallel heights in the ghost layer (Afkhami & Bussmann 2008). The value of C in the ghost layer was also set according to the dynamic contact angle to give a consistent interface normal for the flux computations. In the simulations, we imposed a density ratio of 0.01 to make the simulations stable. We consider this sufficient to reproduce the main features of the water liquid-vapour system.
To evaluate the grid independence of the simulations, we performed a grid refinement study. A refined grid with N x = 512 and N y = 64 was used for θ 0 = 95 • , with Ca = 0.05 and 0.15. The time series of the drop displacements are presented in figure A.1(a). For Ca = 0.05, the drop displacement changed by 1.3% and for Ca = 0.15 by 2.4%. Accordingly, the simulations seem more sensitive to the grid for higher capillary numbers. The difference could appear because of the refinement of the velocity field and because the condition for the dynamic contact angle (2.8) does not completely remove the grid dependence (Legendre & Maglio 2015). Nevertheless, we consider the observed convergence sufficient.
A.2. Customisation for θ 0 = 127 • Using wall-parallel heights in the ghost layer limits the maximum absolute value of the curvature at the contact lines. This curvature corresponds to a maximum forcing on the fluid (|f σ |, equation 2.6). The expression for the curvature is where h = h(y) is the height function giving the x-location of the interface. The sign of κ also depends on the interface orientation; this discussion is not included for brevity (see Aniszewski et al. (2021)). The derivatives are computed with central finite differences. We denote the value of the heights h 0 , h 1 , and h 2 in the ghost, first and second cell layer above the wall, respectively. From the definition of the angle, h 0 = h 1 + ∆ y / tan(θ) (the sign in front of tan(θ) depends on the interface orientation). We assume that h 2 ≈ h 1 . The curvature in the first cell layer becomes 1/ tan(θ) (1 + 1/(4 tan 2 (θ))) 3/2 , (A 2) shown in figure A.1(b) (solid line). The approximation of κ num is zero for θ = 0 • , 90 • , and 180 • . As shown in the figure, |κ num | has one maximum in each of the intervals 0 • < θ < 90 • and 90 • < θ < 180 • . For θ 0 = 127 • , a minimum separation velocity was achieved for a steady-state receding contact line angle smaller than 180 • , as expected from equation (A 2). However, this minimum was not low enough to match the target displacement from MD. To match the MD results, we allowed the usage of the height of the third cell layer, h 3 . The finite-difference scheme of the second derivative was left unchanged. The order of the first derivative, on the other hand, was increased, resulting in 1/ tan(θ) (1 + 1/(9 tan 2 (θ))) 3/2 , (A 3) where we assumed h 3 ≈ h 2 ≈ h 1 . This expression results in significantly higher curvatures for extreme angles (figure A.1b, dashed line). We were then able to match the MD result.
Another possible remedy would be to impose the angle by wall-normal heights. However, in many instances for θ 0 = 127 • , not enough wall-normal heights could be computed at receding nor advancing contact lines for curvature calculations. The current solution is, therefore, more robust. Wall-normal heights could be a viable solution if resolution significantly is increased.

Appendix B. Details of Cahn-Hilliard phase-field model
In this appendix, we provide additional details of the Cahn-Hilliard PF model used in this work. The model is briefly introduced in §2.2. The standard double-well potential is The wetting boundary condition (2.12) contains the so-called switch function, defined as This expression serves as a window function that ensures the wetting boundary condition is applied only at the contact line. Furthermore, the cubic expression (B 2) is not empirical. Instead, it is derived as the equilibrium solution of PF equations based on the selected potential (B 1) and hyperbolic tangent variation of C across the interface. The density and viscosity in the momentum equation (2.2) are spatially dependent. In the PF model, we define these fluid parameters through linear combination based on phase-field variable C as Recall that ρ and µ are the density and the viscosity of the liquid component, and ρ v and µ v are the density and the viscosity of the vapour component.
The introduced governing equations are linearised, written into the weak form, and solved using an open-source finite-element software FreeFEM (Hecht 2012), which allows easy specification of finite-element weak form. Linear elements were used for the phasefield variables, while the fluid flow was resolved using Taylor-Hood elements (quadratic for velocity and linear for pressure). Mesh resolution was refined and results were checked for a few selected simulation cases. The production resolution selected and used for most simulations is ∆s 1 = 3.65 nm far from the interface, and down to ∆s 2 = 0.24 nm within the interface region. The constant time step was used through the simulation as ∆t = 0.002 dimensionless time units. For very small Ca, ∆t was reduced to ensure numerical stability. For simulations with smaller , the mesh was refined at the interface to maintain roughly the same amount of elements over the interface and the time step was reduced to ensure numerical stability. Exact numerical code used to produce the PF results is available freely from Github repository (Lācis & Bagheri 2020-2022. The data to reproduce the figures in the main paper is available in Supplementary File 1.

Appendix C. Details of molecular dynamics simulations
In this appendix, we provide the necessary details for the reader to understand the simulation procedure, interface extraction, and determination of equilibrium angle.

C.1. Numerical implementation
The thin quasi-2D liquid meniscus (figure 1) is composed of 172933 water molecules. The parameters for the SPC/E model are taken from the OPLS-AA force field. The parametrization of SiO 2 quadrupoles is summarized in table C.1. Silicon atoms are treated as virtual sites without mass. Oxygen atoms are restrained to absolute coordinates by a spring of constant κ O . The usage of position restraints grants the substrate some flexibility to re-configure and accommodate water adsorption and desorption. All covalent bonds and angles are treated as rigid constraints. Non-bonded parameters for the interactions between different species are generated via the geometric combination rule. The time-marching step is the same for equilibrium and non-equilibrium runs, δt = 2   fs. We use the leap-frog time marching to update atomistic coordinates. All simulations have been pre-processed and run with GROMACS 2020 (Abraham et al. 2015).
To obtain θ 0 values stated in §3, charge values (q 1 , q 2 , q 3 , q 4 ) = (0.40, 0.60, 0.67, 0.74) e were required. Here, e is the elementary electron charge. Ideal purely-repulsive Lennards-Jones (LJ) walls are placed beyond the silica surfaces at the location of periodic boundary condition. The LJ walls decouple periodic images along y. The starting configuration for production runs is obtained by letting the droplet relax to its equilibrium shape.
The desired shear rate is produced by interpolating position restraints between two reference configurations (figure C.1c). The configuration A is the equilibrium one. In configuration B, horizontal coordinates of the silica layer have been offset by +d ne on the top and by −d ne on the bottom walls. The effective wall velocity is quantified as Here, δλ is the increment of an auxiliary variable λ ∈ [0, ∞) applied at each time step. The λ = 0 corresponds to configuration A and λ = 1 to configuration B. The desired wall velocity is obtained by setting the interpolation increment to δλ = U w δt/d ne .
All simulations are performed in the NVT ensemble at T = 300K and fixing the extent of simulation box to (L x , L y , L z ) = (159.75, 30.634, 4.6765) nm. Periodic boundary conditions are employed along the direction of flow homogeneity z and along the shear direction x, while periodic image interactions along the vertical direction y are avoided by placing ideal Lennard-Jones walls at y = 0 and y = L y (figure C.1a). Bussi-Donadio-Parrinello thermostat (GROMACS 'v-rescale') is applied to both water and silica, with coupling time 0.1 ps for equilibration runs and 10 ps for non-equilibrium runs.
When performing non-equilibrium simulations of liquids one has to bear in mind that most standard choices for thermal coupling either lead to local flow hindering or to artificial cooling where the flow velocity is larger. We estimated the maximum local temperature deviation for the probed range of capillary numbers and we concluded that it has only a marginal effect on steady regimes when surfaces are hydrophilic or moderately hydrophobic. The estimate can be obtained as follows. Imagine a fluid composed of spherical particles. Then temperature and kinetic energy (per particle) can be related at equilibrium via the equipartition theorem where c is the particle's peculiar velocity. For a steady flow and in case of no-slip, molecules close to a solid wall have a deterministic velocity component U w in the x direction. An equilibrium thermostat that is oblivious to hydrodynamics will attempt to re-scale the kinetic energy per particle in order to match the prescribed temperature T 0 , defined as 3 2 The difference between imposed and effective temperature is where Θ is a characteristic temperature differential that tunes how the system is cooled down in the function of the imposed capillary number. For our molecular model, we estimate Θ 0.78 K, which entails that to cool down the near-wall molecules by 1K one needs to prescribe at least Ca 1.132. This rough calculation does not account for the rotational degrees of freedom of water and thus can be regarded as a conservative estimate.
There exist several techniques for correctly thermalizing flow simulations. One among the simplest consists in only coupling the solid substrate and letting the liquid thermalize due to heat exchange. Trying this approach we noticed that in the configurations with θ 0 95 • (which are also problematic due to typically larger Ca) heat transfer between silica and water is not large enough to effectively render the system isothermal. Other techniques would employ either a profile-biased thermostat (Bernardi et al. 2010) or a dissipative-particle-dynamics thermostat (Soddemann et al. 2003;Goga et al. 2012). However, these thermostats are not currently implemented in GROMACS.
Hydrodynamic fields (density, velocity, and temperature) are directly measured from MD trajectories. Each quantity is averaged in space on a grid with spacing (h x , h y ) (0.20, 0.20) nm (figure 1b), and over time by aggregating all measurements in consecutive windows of 12.5 ps. Averaged and binned variables are saved to file "on-the-fly" concurrently with the simulation, thus vastly reducing the output size. Saving all atomistic trajectories would not be feasible. Consequently, the output of MD simulations is a range of data files containing so-called frames, corresponding to the sequence in time of the partially averaged MD data. Each frame contains the instantaneous field outputs ρ i (x, y, t), u i x (x, y, t) and u i y (x, y, t) as defined in §3. Post-processing to obtain averages over time intervals of several ns has been performed with in-house codes based on the freely available repository (https://github.com/MicPellegrino/densmap.git). The exact scripts are available upon reasonable request.

C.2. Equilibration runs and centre of mass correction
To determine that the signature of the initialisation is fully disappeared, we first visually inspect the time series of the relaxing contact angle. Based on the decay of the signal we have determined conservative cut-off times for different q values. Then we check a-posteriori the cross-correlation between the signals at each contact line. This is done to ensure that any transient relaxation dynamic has disappeared and that the size of the molecular system is large enough to effectively localize contact line motion. After the cut-off time, we have continued the runs to collect a sufficient amount of statistics for the measurements, at least 4 ns for all q values.
From MD simulations of the equilibrium configuration, we obtain many sequential frames of hydrodynamic variables. To measure the geometrical features, such as the local interface curvature, we average all frames in the equilibrium state by shifting the centre of mass (COM) of the liquid droplet in each frame to the centre of the domain. The reason behind this procedure lies in the fact that COM correction is turned off in the MD simulations themselves. Run-time COM correction can potentially hinder relaxation in equilibrium simulations and create velocity measurement artefacts in nonequilibrium ones. This averaging procedure is employed before interface extraction for both equilibration and sheared MD runs.

C.3. Interface extraction and θ 0 measurement
In this appendix, we describe the extraction of the interface shape from the water density distribution ρ(x, y). The distribution for equilibration run with θ 0 = 95 • along the bin with vertical coordinate y ≈ 1.5 nm is shown in figure 5(e,f). We consider two criteria to define the exact x ni coordinate of the interface, i.e.
(C 5a,b) These are based on global liquid density (C 5a) and slice liquid density (C 5b) at an n-th vertical bin with coordinate y n . If the sought density value is not located in any single bin, linear interpolation is used to find the exact coordinate. The interface extraction according to (C 5a) we call "global interface extraction" and the obtained interface we call "global interface". This is the approach typically used in literature. However, the density layering (figure C.2b) exhibits itself in the interface shape near the surface (figure C.2a, green line). The interface extracted according to (C 5b) we call "slice interface". The x ni obtained from (C 5b) show reduced layering influence (figure C.2a, blue line). Further away from the wall results from (C 5a,b) agree. From (C 5b) it is also possible to define the interface point for adsorbed water layer (figure C.2a, y = 0.7 nm). However, this is only a visual evaluation of both extraction procedures and it is not yet clear if any one of those is advantageous when comparing MD with CFM. As the next step, we measure the equilibrium contact angle θ 0 . We use both interface shapes extracted according to (C 5). From the interface shape, we compute the interface angle θ(y) along the height of the interface, as defined in figure 2(a). The angle is obtained from the slope of the interface segments (encircled with green in figure C.2a). Therefore angle measurements are located at the boundaries of the MD bins. The Young-Laplace equation for constant surface tension is where ∆p is the pressure jump across the interface and γ is the local curvature of the interface. In equilibrium, we expect the pressure in the whole droplet to be constant and the pressure in the vapour phase to be negligible. According to the Young-Laplace where c 1 and c 2 are constants determined by boundary conditions and s is curvilinear coordinate along the interface. For convenient comparison, we transfer θ(y) to θ(s). Since the theoretical function is known, we fit the obtained MD results with (C 7). We observe that both global and slice interfaces rapidly deviate from linear relationships near the wall (figure C.2d). The same effect is observed in figure C.2(a), where the interface segments closest to the wall exhibit significantly different angles compared to segments ≈ 1.5 nm above the wall. Therefore we conclude that several interface points closest to the wall can not be used for a reliable comparison with continuum description. Next, we determine how many MD points near the wall need to be excluded. We do this by gradually removing the closest interface points near the wall. The procedure is carried out until the standard error of the linear fit (C 7) reaches its minimum. This process is applied to both global and slice interfaces. We mostly observe that a larger number of points have to be removed from the global interface shape to obtain the minimum in the standard error. Therefore for producing MD results in this work, we choose to focus only on the slice interface. We postulate that the first remaining point on the interface is the first reliable interface measurement for comparison with CFM. The obtained θ (s) from slice interface for the remaining interface points is presented in figure C.2(c) with a blue line. The corresponding best linear fit is given with a red dotted line. To obtain the equilibrium angle θ 0 , we use the hydrodynamic wall position assumed in the main paper. The wall position is shown in red in figure C.2(b). We extrapolate the linear fit to the assumed positions. The equilibrium angle is computed as an average of extrapolated values at the top and bottom walls (figure C.2c, red crosses). However, the total arc length of the interface depends on the contact angle. Therefore the equilibrium angle θ 0 is iteratively obtained in the following procedure. We first centre the MD data to ideal arc length for given θ 0 . Then the agreement with the given θ 0 is verified by extrapolating the best linear fit of the MD data to wall locations. Finally, the next θ 0 is set as the averaged of the previous estimate and the current estimate. The final obtained equilibrium angle in figure C.2(c) is presented with black dotted lines. Inset (figure C.2d) shows the angle deviation magnitude from the linear expression near the wall. The deviation is much larger than the noise observed in the bulk of the MD results. Consequently, the deviation near the wall can not be explained by the thermal fluctuations of MD. Therefore some other molecular effects are in play.
The process of removing the unreliable MD interface points and determining the equilibrium contact angle θ 0 is repeated for all considered MD surface charges. Obtained first reliable interface extraction bin locations for comparison with CFM are summarized in the main paper figure 5(a-d) and the measured contact angles θ 0 are 127 • , 95 • , 69 • and 38 • as stated in the main paper text.

C.4. Polynomial extrapolation of MD interface shape
In this appendix, we describe how a polynomial fit is used to determine the final drop displacement for comparison with CFM models. In addition, the usefulness of the fit for extracting dynamic contact angle is assessed. Recall that through the extraction of θ 0 values (appendix C.3) the first reliable bins for comparison with CFM have been identified (shown with green in figure 5a-d). Consequently, there is a gap in interface data near the walls. To illustrate this, we show the equilibrium and sheared (Ca = 0.20) interface shapes for θ 0 = 95 • in figure C.3(a). With a red dashed line, we show the assumed wall position. In CFM models, however, the interface shape continues to evolve smoothly until meeting the wall. Consequently, when comparing MD with CFM models, the empty data region near the wall is undesirable.
To remedy this issue, we introduce a fit of the interface shape. Unlike for the equilibration run -for which the function corresponding to the interface shape was known -, the functional form of the dynamic interface is not known. Therefore we choose a polynomial with some order N p that is not yet known. To fix the N p parameter, a convergence study of extrapolated drop displacement by varying the N p parameter is carried out. The example result of the convergence study is shown in figure C.3(c) for θ 0 = 95 • and Ca = 0.20. We typically observe that initially there are large changes of the extrapolated drop displacement. However, after some order, the magnitude of difference reduces. To settle on the final polynomial order, we use the following guidelines. As a rule of thumb, we choose the order after which the drop displacement ∆x visually seems to oscillate around some value. In addition, we introduce an arbitrary limit N p < 12 to avoid over-fitting the MD data. Finally, we evaluate the agreement with the interface angle (figure C.3b) and increase the order if the agreement is not satisfactory. The final chose polynomial order for Ca = 0.20 is N p = 7. The final obtained steady displacement from MD is ∆x = 13.89 nm. The fitted polynomial for the interface, the shape is shown in figure C.3(a,b) with a magenta dotted line.
While the ∆x can be obtained by extrapolating the polynomial fit, this approach does not give reliable measurements of the dynamic contact angle at the wall. To illustrate this, in figure C.3(b) we add two more polynomial fits with N p = 8 and N p = 11. All three polynomial orders give very similar ∆x values (C.3c). However, the predicted advancing dynamic contact angles at the wall (magenta crosses in figure C.3b) are significantly different. Similar differences have also been observed for receding contact lines for other simulations. Consequently, the first reliable measurement of the dynamic contact angle can be taken only at some distance from the wall (green crosses in figure C.3b). The polynomial extrapolation in the main paper is used to get a more accurate a posteori steady ∆x measurement and to read off the dynamic contact angle at the reliable bin location. The convergence of ∆x is rechecked for each unsteady MD simulation by following the guideline explained above.

Appendix D. Choice of vapour properties
According to empirical formulae (Engineering ToolBox 2004), the water vapour saturation pressure (which is the same as the gas pressure in the absence of other gasses) at T = 300 K is The water vapour density is The viscosity of water vapour (also called steam) can be looked up in tables (Engineering ToolBox 2014). Linear interpolation between two given values closest to T = 300 K gives us These are the parameters reported in the main paper, figure 1. It was also checked that the results are only weakly sensitive to the exact value of the vapour viscosity.

Appendix E. Wall location and no-slip condition
In this appendix, we motivate the choice of the hydrodynamic wall position and the applicability of the no-slip condition. We also how that small changes in wall location can influence results near the contact line significantly. Recall that the wall position in CFM is set at the centre of the bin with coordinate y = 0.7 nm (black line in figure E.1a). If the chosen wall location is appropriate, the CFM should predict flow velocity accurately even a very small distance above the wall. Therefore we select the next bin with coordinate y = 0.9 nm (green bin in figure E.1a) to compare velocity distribution between the MD and PF models. We omit VOF from the comparison for clarity and we look at all calibration simulations. We extract u x from MD at the green bin. The MD velocity data is obtained with help of two averaging approaches. The first approach is the MD frame average over the steady regime together with correction for the COM, the same way as done in appendix C.2. This provides the global flow field data and interface shape. Locally, as we approach the two-phase interface, the stochastic interface oscillations become present and can influence the measurement of the hydrodynamic variables. To reduce this influence, we repeat the averaging procedure over all frames in a steady regime, but instead of centring those around the COM, we centre them around the instantaneous interface positions at the left and the right side of the drop. Using this approach, we obtain a cleaner signal from MD closer to the interface. To obtain a single velocity curve, we move the interface averaged results to the global interface location. Then, we replace the velocity data from the COM averaged data with the interface average data until ≈ 10 nm away from the contact line. The noisy data on the vapour side is neglected. Finally, to further reduce the noise in MD, we make use of the symmetry in the system and take the mean between profiles obtained at the bottom and top walls. The obtained MD stream-wise velocity distribution is shown in figure E.1(b-e) with solid black lines. Contact angles and Ca numbers are presented in the title of individual panels.
For comparison, we extract the u x distribution along x coordinate at the same y location from all calibration no-slip PF simulations. The PF results for θ 0 = 127 • are shown in figure E.1(a) with a green solid line. We observe that PF results do not agree with MD results over the full span of x coordinate. We repeat the PF simulation, by gradually increasing the s value until a good match is obtained. Through this, we obtain s = 0.44 nm. The PF velocity results with s = 0.44 nm are shown in figure E.1(a) with a solid red line. Now, a good agreement between MD and PF is obtained.
For θ 0 95 • (figure E.1c-e) we observe that agreement between PF velocity predictions and MD results is good at s = 0. Very good correspondence is obtained at the centre of the drop. As contact line regions are approached, the agreement deteriorates. However, the agreement below the liquid bulk is sufficient to conclude that the no-slip condition ( s = 0) is appropriate for these contact angles. The hydrodynamic wall position in the current work is essentially an assumption. Therefore we have also investigated the sensitivity of the velocity profile obtained from PF to small perturbations in wall location. With grey area (figure E.1b-e) we show the region in which the PF results would fall if the wall location would be moved up or down by 0.1 nm. These results are obtained by sampling by 0.1 nm closer or further away from the solid wall. It was verified that this is equivalent to actually changing the wall position and also the channel height. By looking at the grey region, we observe that the variations are very narrow for all θ 0 values in the centre of the drop. However, approaching the contact line region for θ 0 = 95 • and 69 • the variation grows. This suggests that for the description of the processes near the contact line the exact location of the solid wall could play a significant role.
It has to be recognised that the slip length of MD systems has been studied extensively before, see for example works by Thompson & Troian (1997) and Huang et al. (2008). In particular, Huang et al. (2008) investigated the slip length of SPC/E water over a range of surfaces, from silane monolayers to more common Lennards-Jones models. They found that, up to a good accuracy, s ∼ (1 + cos θ 0 ) −2 . It has to be marked that, to the best of authors knowledge, similar study in MD with surfaces that form hydrogen bonds with water has not been carried out before and is out of scope of the current study. Nevertheless, the slip length values obtained in current work qualitatively agree with results of Huang et al. (2008) -as θ 0 increases, the slip length grows as well. However, the current inaccuracy of wall location and selected binning resolution (0.2 nm) prohibits determination of the exact slip length variation for contact angles θ 0 = 38 • − 95 • .
Finally, it is interesting to note that the liquid density variations ( §3) do not impede reliable flow velocity measurements. Reliable velocity measurements can be taken closer to the wall if compared to reliable interface angle measurements (compare red and green bin locations in figure 5a-d). In addition, through determining the validity of the noslip condition, we have demonstrated that PF can accurately predict the liquid velocity distribution as close as 0.2 nm above the last oxygen atom of the solid substrate.

Appendix F. Streamlines from PF near the contact line
To deepen the understanding of the parameter M , µ f and influence on the PF results, the steady flow field in the vicinity of the receding contact line is investigated. For this, the PF simulation with θ 0 = 95 • and Ca = 0.20 is run until the steady ∆x is reached. The flow field at the last time instant is used to compute the streamlines. With black lines in figure F.1 we present streamlines near the bottom left receding contact line in a zoomed-in window of roughly 12 nm × 7 nm. The two-phase interface, defined as C = 0, is presented with a thick red line. The thick blue line identifies a streamline that originates from within the liquid drop at a 0.5 nm distance from the bottom wall. This particular streamline can be leveraged to compare the amount of streamline crossing and the amount of overshoot at the two-phase interface.
In the first row (figure F.1a-c) the influence of PF mobility is shown. First, we investigate the flow field with the smallest M ( figure F.1a). Overall, the streamlines are similar to the wedge solution derived and presented by Huh & Scriven (1971) and the PF solution thoroughly analysed by Seppecher (1996). On the large viscosity side (in the liquid part), there is only one vortex, while on the small viscosity side (in the vapour part) there are two adjacent vortices. Due to the diffuse nature of the model, the stagnation point is displaced slightly to the left and above the contact line at the wall defined by C = 0. In addition, it can be observed that the blue streamline approaches the two-phase interface, then turns and follows the two-phase interface tangentially. By setting M ten times larger (figure F.1b), it can be observed that the streamline crossing over the interface is increased. The blue streamline now crosses the interface (red line) and continues in the vapour phase. Whereas streamline originating at a slightly larger distance from the wall (see a black line that partially overlaps with the blue line) turns and follows the interface tangentially in the vapour side. For M hundred times larger (figure F.1c), the streamline crossing is increased even more. The blue streamline proceeds straight into the vapour phase, and so do the streamlines originating up to a distance of roughly 2 nm above the wall. The original wedge flow pattern can barely be recognised.
The second row (figure F.1d-f) illustrates the effect of varying contact line friction µ f . When there is no friction (figure F.1d), the contact angle at the wall is equal to θ 0 and the marked streamline overshoots the interface by around 1 nm and is pulled back within the liquid drop higher above the wall. Increasing friction to µ (figure F.1e) leads to an overshoot of the streamline around 2 nm. Higher above the wall, the streamline is pulled back and continues parallel to the interface at roughly the same distance from the interface as observed in figure F.1(d). The contact angle at the wall deviates from θ 0 . For the largest contact line friction (figure F.1f) we observe an even more pronounced contact angle departure from θ 0 . In addition, the marked streamline crosses the interface and continues in the vapour phase. Finally, in the third row (figure F.1g-i) we exemplify the effect of the interface thickness. As is reduced, the behaviour of the blue streamline changes from full crossover to vapour phase ( figure F.1i), to a small overshoot of an order of nm (figure F.1h) and finally to no crossing over the interface (figure F.1g). The interface shape, on the other hand, remains practically the same for = 0.70 nm and 0.18 nm, which is again a signature of the sharp interface limit.