1. 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 ends up naturally with an immobile line at which both fluid phases meet with the solid, the 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 (Reference Huh and Scriven1971) as a stress singularity at the contact line. Motivated by the importance of the moving contact line in applications such as printing (Kumar Reference Kumar2015), 
${\rm CO}_2$ storage and water management in fuel cells (Singh et al. Reference Singh, Jung, Brinkmann and Seemann2019), 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. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013; Sui, Ding & Spelt Reference Sui, Ding and Spelt2014) 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 Reference Voinov1976; Cox Reference Cox1986; Shikhmurzaev Reference Shikhmurzaev1994, Reference Shikhmurzaev1997; Kalliadasis, Bielarz & Homsy Reference Kalliadasis, Bielarz and Homsy2000; Pismen & Pomeau Reference Pismen and Pomeau2000; Eggers Reference Eggers2004; Flitton & King Reference Flitton and King2004; Snoeijer Reference Snoeijer2006; Wilson et al. Reference Wilson, Summers, Shikhmurzaev, Clarke and Blake2006; Pismen & Eggers Reference Pismen and Eggers2008; Snoeijer & Eggers Reference Snoeijer and Eggers2010; Chan, Snoeijer & Eggers Reference Chan, Snoeijer and Eggers2012; Nold et al. Reference Nold, González MacDowell, Sibley, Goddard and Kalliadasis2018; Chan et al. Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020), a theoretical consensus is yet to be reached. To a large extent, progress is hindered by a lack of understanding of the nanoscale physics of wetting phenomena (Afkhami, Gambaryan-Roisman & Pismen Reference Afkhami, Gambaryan-Roisman and Pismen2020; Afkhami Reference Afkhami2022). For example, a fundamental question relates to the nanoscopic contact angle. The typical choice until now has been a constant equilibrium angle (Kronbichler & Kreiss Reference Kronbichler and Kreiss2017). However, recent experimental (Deng et al. Reference Deng, Chen, Liu, Yu and Wang2016) and numerical (Fernández-Toledano, Blake & De Coninck Reference Fernández-Toledano, Blake and De Coninck2021) evidence suggests a contact line velocity-dependent dynamic contact angle even at nanoscale. The development of more accurate measurement techniques is ongoing (Thormann et al. Reference Thormann, Simonsen, Hansen and Mouritsen2008; Eriksson et al. Reference Eriksson2019); 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, Ogonowski & De Coninck Reference Gentner, Ogonowski and De Coninck2003; Smith et al. Reference Smith, Müller, Craster and Matar2016, Reference Smith, Theodorakis, Craster and Matar2018; Perumanath et al. Reference Perumanath, Borg, Chubynsky, Sprittles and Reese2019; Lācis et al. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020), 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 model numerically 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 Reference Navier1823; Spelt Reference Spelt2005) or by an implicit numerical slip at the contact line (Renardy, Renardy & Li Reference Renardy, Renardy and Li2001; Afkhami, Zaleski & Bussmann Reference Afkhami, Zaleski and Bussmann2009). For diffuse interface models, such as the Cahn–Hilliard phase-field (PF) model (Jacqmin Reference Jacqmin2000), 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 (Reference Shan and Chen1993) or Gunstensen et al. (Reference Gunstensen, Rothman, Zaleski and Zanetti1991) models. It is argued that in the Shan–Chen model, the contact line moves through phase change (Kamali et al. Reference Kamali, Gillissen, Sundaresan and Van den Akker2011), similar to the LBM with thermodynamically consistent potentials (Briant, Wagner & Yeomans Reference Briant, Wagner and Yeomans2004). The LBM has also been used to solve PF equations (Briant & Yeomans Reference Briant and Yeomans2004); 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 (Reference Legendre and Maglio2015). 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 Reference Yue and Feng2011). A part of the calibration is choosing the diffuse interface thickness in a manner to satisfy the so-called sharp interface limit (Yue, Zhou & Feng Reference Yue, Zhou and Feng2010; Magaletti et al. Reference Magaletti, Picano, Chinappi, Marino and Casciola2013; Xu, Di & Yu Reference Xu, Di and Yu2018). 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, Wang & Sheng Reference Qian, Wang and Sheng2003; Barclay & Lukes Reference Barclay and Lukes2019; Mohand et al. Reference Mohand, Hoang, Galliero and Legendre2019). The typical MD work considers Lennard-Jones types of surfaces with large slip lengths. Only more recently have water MD simulations been carried out by some of us (Johansson & Hess Reference Johansson and Hess2018; Lācis et al. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020) over surfaces with negligible slip.
In this work, we generate benchmark data from MD simulations of wetting over no-slip surfaces. These conditions can be reproduced employing simple-point-charge/extended (SPC/E) water on a smooth silica-like substrate (Johansson, Carlson & Hess Reference Johansson, Carlson and Hess2015). We choose a forced wetting set-up (Blake et al. Reference Blake, Fernandez-Toledano, Doyen and De Coninck2015) – instead of a capillary-driven one (Villanueva & Amberg Reference Villanueva and Amberg2006) – 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 Reference Jacqmin2004; Sbragaglia, Sugiyama & Biferale Reference Sbragaglia, Sugiyama and Biferale2008; Gao & Lu Reference Gao and Lu2013; Wang & McCarthy Reference Wang and McCarthy2013) 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 droplet configuration provides rich interfacial dynamics that are challenging to capture with CFM models.
(a) Dimensions and properties of the sheared-droplet configuration considered in the current work. (b) Close-up of the molecular system near the moving contact lines. Overlying solid mesh illustrates the binning boundaries for the collection of flow data from MD simulations.
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 (
$\theta _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,\theta _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 an extension of the work presented by Lācis et al. (Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020), 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 MD of the sheared droplet reported herein provide insight into 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 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 and 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–F, important physical and technical details are provided.
2. 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 horizontal plates located at 
$y=0$ and 
$y=H$. A liquid drop of density 
$\rho _\ell$ and viscosity 
$\mu _\ell$ is sandwiched between the plates such that its maximum width is 
$W$. The drop is surrounded by water vapour with density 
$\rho _v$ and viscosity 
$\mu _v$. The surface tension between the phases is constant and denoted by 
$\sigma$. 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 
$\theta _0$. The former is discussed interchangeably in its non-dimensional form, using the capillary number
\begin{equation} {Ca} = \frac{2\mu_\ell U_w}{\sigma}. \end{equation}
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, 
$\theta _0$, is a measure of the surface's affinity to water. In hydrophilic conditions (
$\theta _0<90^\circ$), the affinity is strong, resulting in a larger droplet deformation compared to hydrophobic conditions (
$\theta _0>90^\circ$). Under dynamic conditions, the contact line can be different from the equilibrium one. For advancing contact lines (marked ‘A’ in figure 2a), the liquid displaces the vapour, while for receding contact lines (marked ‘R’ in figure 2a), the vapour displaces the liquid. Both these processes are largely determined by molecular interactions between the water and the substrate as depicted in figure 1(b).
(a) Drop displacement measurement at the left (
$\Delta x_l$) and right (
$\Delta x_r$) sides at selected time instances. The interface angle is defined with respect to the horizontal line and is sampled as a function of the vertical coordinate 
$\theta (y)$. Note that 
$\theta (y)$ is measured only for 
${Ca} < {Ca}_{c}$, while 
$\Delta x$ is measured for all 
$U_w$. (b) Variations of drop displacement 
$\Delta x = ( \Delta x_l + \Delta x_r )/2$ over time in VOF simulations for 
$\lambda = 0.47, 0.66, 0.94, 1.87, 3.74$ nm (increasing with the arrow). The equilibrium contact angle is 
$\theta _0 = 95^\circ$, and the capillary number is 
${Ca} = 0.20 < {Ca}_{c}$.
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 the left (
$\Delta x_l$) and right (
$\Delta x_r$) two-phase interface, respectively. For more detailed characterization, we also evaluate the interface angle 
$\theta (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 
$(\boldsymbol {{u}}, p)$ are obtained by solving the incompressible Navier–Stokes equations in the domain containing both phases. In two dimensions, the equations are
$$\begin{gather} \rho ( x,y ) \left[ \frac{\partial \boldsymbol{{u}}}{\partial t} + \left(\boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{{u}} \right] ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\mu ( x,y ) \left\{ \boldsymbol{\nabla} \boldsymbol{{u}} + (\boldsymbol{\nabla} \boldsymbol{{u}} )^{\rm T} \right\}\right] + \boldsymbol{{f}}_\sigma, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{u}} = 0. \end{gather}$$
Here, 
$\boldsymbol {{f}}_\sigma$ is the surface tension force, and 
$\rho ( x,y )$ and 
$\mu ( x,y )$ are spatially dependent density and viscosity, respectively. The functions 
$\rho ( x,y )$ and 
$\mu ( 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 
$\boldsymbol {{f}}_\sigma$ 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
\begin{equation} u_x = u_w + \ell_s\,\frac{\partial u_x}{\partial n}, \end{equation}
where 
$\ell _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 
$\ell _s = 0$ for the chosen liquid–solid combination, the implementation of the CFM models allows for tests with non-zero value.
Equations (2.2) and (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 the two-phase interface are treated differs.
2.1. 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
\begin{equation} \frac{\partial C}{\partial t} ={-} \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C. \end{equation}The surface tension force is applied by the continuous surface force (CSF) method (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992):
\begin{equation} \boldsymbol{{f}}_\sigma ={-}\sigma \kappa \boldsymbol{\nabla} C, \end{equation}
where 
$\kappa$ is the curvature of the interface.
As a boundary condition for the phase variable, we impose a dynamic contact angle 
$\theta _{num}$. The contact angle is the angle between the interface (
$\hat {n}_i$) and wall (
$\hat {n}$) normals, so
\begin{equation} \hat{n}_i\boldsymbol{\cdot} \hat{n} = \cos(\theta_{num}), \quad \text{where} \quad \hat{n}_i ={-}\frac{\boldsymbol{\nabla} C}{|\boldsymbol{\nabla} C|}, \end{equation}
and the sign might vary depending on convention. The relationship for 
$\theta _{num}$ is inspired by Cox's theory, as described and evaluated by Legendre & Maglio (Reference Legendre and Maglio2015). It takes the form
\begin{equation} G(\theta_{num}) = G(\theta_0) + {Ca}_{cl}\ln\left(\frac{\varDelta/2}{\lambda} \right), \end{equation}
where 
$G$ is a monotonically increasing function, 
$\lambda$ is a microscopic cut-off length scale, and 
$\varDelta$ is the wall-normal cell height. We confirm the grid convergence reported by Legendre & Maglio (Reference Legendre and Maglio2015) 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 interpolating the velocity field linearly at the first grid cell. For the most hydrophobic and hydrophilic configurations, only the constant wall velocity was used for increased robustness.
The geometric VOF model is solved with an open-source finite volume code, called the PArallel, Robust, Interface Simulator (PARIS) (Aniszewski et al. Reference Aniszewski2021). The full set of equations (2.2), (2.3), (2.5) and (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 (
$\theta _0 = 127^\circ$), 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 
$\lambda$. To test the effect of 
$\lambda$, we select 
$\theta _0 = 95^\circ$ and 
${Ca} = 0.20 < {Ca}_{c}$. The evolution of the drop displacement 
$\Delta x$ in time for 
$\lambda$ 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, 
$\Delta x_l = \Delta x_r = \Delta x$ in the CFM simulations, which we have verified numerically. We observe that larger 
$\lambda$ correspond to smaller steady 
$\Delta x$. As 
$\lambda$ is increased, the contact angle at the surface (2.8) and thus the interface shape (2.7) near the wall are 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 
$\Delta x$. The 
$\Delta x$ for the initial time agrees between all 
$\lambda$ values. In this short period, the contact line does not slip with respect to the substrate, and the slope of 
$\Delta x (t)$ is equivalent to the wall separation velocity (
$2 U_w$).
2.2. Cahn–Hilliard PF model
We choose a model of a binary mixture with a classical fourth-order polynomial potential 
$\varPsi$ (see (B1) in Appendix B) (Jacqmin Reference Jacqmin2000; Carlson Reference Carlson2012). 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
\begin{equation} \frac{\partial C}{\partial t} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ M \boldsymbol{\nabla} \phi \right] - \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C . \end{equation}
Here, 
$M$ is the PF mobility and 
$\phi$ is the chemical potential. The latter is defined as
\begin{equation} \phi = \frac{2\sqrt{2}}{3}\,\frac{\sigma}{\epsilon}\,\varPsi'(C ) - \frac{2\sqrt{2}}{3}\,\sigma\epsilon \nabla^2 C. \end{equation}
The chemical potential (2.10) contains the surface tension (
$\sigma$) and the interface thickness (
$\epsilon$). The derivation of (2.9) and (2.10) is standard and can be found in Carlson (Reference Carlson2012) and Jacqmin (Reference Jacqmin2000). The surface tension 
$\sigma$ is a physical parameter and is set to a value corresponding to the water liquid–vapour interface. The constants 
$M$ and 
$\epsilon$, on other hand, are typically treated as numerical parameters. Having solved for function 
$C$, the surface tension force in (2.2) is inserted as
\begin{equation} \boldsymbol{{f}}_\sigma = \phi\, \boldsymbol{\nabla} C. \end{equation}The Cahn–Hilliard equation (2.9) is a fourth-order partial differential equation, thus two boundary conditions are needed on solid walls. The lowest-order boundary condition is
\begin{equation} - \mu_f \epsilon \left( \frac{\partial C}{\partial t} + \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C \right) = \frac{2\sqrt{2}}{3}\,\sigma\epsilon\,\boldsymbol{\nabla} C \boldsymbol{\cdot} \hat{n} - \sigma \cos ( \theta_0 )\,g'( C ), \end{equation}
where 
$\mu _f$ is contact line friction, and 
$g$ is switch function (B2). Equation (2.12) is also known as the non-equilibrium wetting condition and requires the equilibrium contact angle 
$\theta _0$ as an input. Setting a non-zero 
$\mu _f$ yields a dynamic contact angle different from 
$\theta _0$ (Jacqmin Reference Jacqmin2000; Qian et al. Reference Qian, Wang and Sheng2003; Carlson, Do-Quang & Amberg Reference Carlson, Do-Quang and Amberg2011). The second boundary condition on solid walls is 
$\boldsymbol{\nabla} \phi \boldsymbol {\cdot } \hat {n} = 0$, which states that there is no diffusive flux through the walls.
The full system of fluid and PF equations (2.2), (2.3), (2.9) and (2.10), together with the boundary conditions, is discretized and solved using open-source code FreeFEM (Hecht Reference Hecht2012). 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 PF mobility 
$M$, the contact line friction 
$\mu _f$, and the interface thickness 
$\epsilon$. To provide an intuition about each parameter, we vary them one at a time. As before, 
$\theta _0 = 95^\circ$ and 
${Ca} = 0.20$. Figure 3(a) shows 
$\Delta x$ for 
$\mu _f = 0$, 
$\epsilon = 0.7$ nm and 
$M = (3.5; 7.0; 10.5; 14.0; 17.5) \times 10^{-16}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. We observe that the final steady 
$\Delta x$ value is reduced as 
$M$ is increased, i.e. smaller deformation with increased diffusion. In figure 3(b), the displacement for 
$\epsilon = 0.7$ nm, 
$M = 1.75 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and 
$\mu _f = (0.0; 0.5; 1.0; 1.5; 2.0)\,\mu _\ell$ is shown. Here, 
$\Delta x$ increases with 
$\mu _f$, therefore the contact line friction has an opposite effect compared to mobility. This competition has been previously clearly showcased by Yue & Feng (Reference Yue and Feng2011). Finally, the evolution of 
$\Delta x$ in time for 
$\epsilon = (0.18;0.35;0.70;1.40;2.80 )$ nm is shown in figure 3(c) with 
$\mu _f = 0$ and 
$M = 1.08 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. It can be observed that – in contrast to the 
$M$ and 
$\mu _f$ variations – notable differences are present for initial times. As 
$\epsilon$ is decreased (in the direction against the arrow), the steady 
$\Delta x$ value is converging. This is a signature of the sharp interface limit (Yue et al. Reference Yue, Zhou and Feng2010; Xu et al. Reference Xu, Di and Yu2018). For enriched understanding, in Appendix F we report streamlines from PF near a receding contact line for similar parameter variations.
Drop displacement 
$\Delta x$ over time in PF simulations, varying PF mobility (a), contact line friction (b), and interface thickness (c). In (a), 
$\mu _f = 0$, 
$\epsilon = 0.7$ nm and PF mobility is 
$M = (3.5; 7.0; 10.5; 14.0; 17.5) \times 10^{-16}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. In (b), 
$\epsilon = 0.7$ nm, 
$M = 1.75 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and contact line friction is 
$\mu _f = (0.0; 0.5; 1.0; 1.5; 2.0)\,\mu _\ell$. In (c), 
$\mu _f = 0$, 
$M = 1.08 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and interface thickness is 
$\epsilon = (0.18;0.35;0.70;1.40;2.80 )$ nm. Varying parameters are increasing along the black arrows in all panels. In all simulations, the equilibrium contact angle is 
$\theta _0 = 95^\circ$ and the wall velocity is 
$U_w = 6.67\ {\rm m}\ {\rm s}^{-1}$.
3. 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 
${\rm SiO}_2$ quadrupoles (figure 4a) that are restrained into a hexagonal lattice. A quasi-two-dimensional 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 Reference Johansson and Hess2018). 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.
(a) Molecular geometry of silica quadrupoles. (b) Water molecules adsorbing (bonds sketched with dashed lines) and forming hydrogen bond (sketched with dotted line) with silica substrate. (c) Top view of the contact line region, showing the lattice of silica quadrupoles and water up to 
$\approx$1.5 nm above the lower periodic boundary (blue line in b).
The strength of the water–substrate interaction can be tuned by adjusting the charge distribution in 
${\rm SiO}_2$ (figure 4a). Different interaction strengths will result in different equilibrium contact angles 
$\theta _0$. We simulate the system via atomistic MD in the NVT ensemble, using well-established force fields and thermostats. Details regarding the physical and numerical simulation set-up can be found in Appendix C.1.
3.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 so-called equilibration runs (see Appendix C.2). During the run time of the MD simulation, we collect flow data in regular 
$0.2\ {\rm nm} \times 0.2\ {\rm nm}$ bins for a time interval 12.5 ps. This yields instantaneous density 
$\rho ^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, 
$\rho ^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. 
$\rho (x,y)$. The interface position is determined by seeking the location where the liquid density transitions from bulk density to very small vapour density. Examples of density distributions are shown in figures 5(e, f). More details on interface extraction are provided in C.3.
Liquid water density variation near the bottom wall for equilibrium angles (a) 
$\theta _0 = 127^\circ$, (b) 
$95^\circ$, (c) 
$69^\circ$, and (d) 
$38^\circ$. Dashed black lines illustrate the boundaries of the bins. The bin filled with red shows the selected location of the solid wall, while the bin filled with green shows the first reliable interface measurement. Water density distributions along the green bin for 
$\theta _0 = 95^\circ$ configuration are shown (e) over the full span of 
$x$ coordinates, and (f) near the left liquid–vapour interface.
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 in figures 5(a–d). The 
$q$ values are tuned to yield 
$\theta _0 = 127^\circ$, 
$95^\circ$, 
$69^\circ$ and 
$38^\circ$. This allows us to investigate hydrophobic, neutral and hydrophilic wetting conditions. In parallel to 
$\theta _0$ extraction, we also identify the first reliable bin for interface shape measurement. This bin is shown in green in figures 5(a–d). The interface angle computed from points closer to the wall exhibits extreme deviation from continuum description (figure 20d). 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 (L-J) liquids near solid walls and two-phase interfaces (Bugel, Galliéro & Caltagirone Reference Bugel, Galliéro and Caltagirone2011; Stephan et al. Reference Stephan, Liu, Langenbach, Chapman and Hasse2018). We determine the extent of 
$\rho$ oscillations near the wall in our MD system. The water liquid density distribution along the height of the channel is computed as
\begin{equation} \rho_y (y) = \int\limits_{x_l}^{x_r} \rho (x,y)\,{{\rm d}\kern0.06em x}. \end{equation}
Here, the boundaries for integration 
$x_l$ and 
$x_r$ are selected for each 
$\theta _0$ to fall within the liquid phase for all 
$y$ coordinates. The close-up of 
$\rho _y (y)$ near the bottom wall is shown in figures 5(a–d). Oscillations in 
$\rho$ have smaller amplitude and occur over smaller distances than observed typically in L-J systems. The small layering is an outcome of the combination of the SPC/E water model and the 
${\rm SiO}_2$ surface model. For the same value of 
$\theta _0$, the layering would propagate further into the liquid phase if an L-J surface were used instead of 
${\rm SiO}_2$.
3.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 
$\Delta x$ is obtained as an average between the left (
$\Delta x_l$) and right (
$\Delta x_r$) interfaces. We determine the drop displacement using the first reliable bins near the top and bottom walls (figure 5, green bins). Since there are 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 
$\Delta x (t)$ for 
$\theta _0 = 95^\circ$ for different capillary numbers, starting at 
${Ca} = 0.05$ and then increased incrementally by 
$\Delta {Ca} = 0.05$. For all 
${Ca}$ up to and including 
${Ca} = 0.25$, we observed a stable configuration. The obtained steady drop displacements for 
$\theta _0 = 95^\circ$ are summarized in the 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 three simulations in a stable regime are gathered. For all 
$\theta _0$, as expected, we observe that 
$\Delta x$ increases with 
${Ca}$.
Time evolution of 
$\Delta x$ from MD with 
$\theta _0 = 95^\circ$ for simulations reaching steady state. In (a), 
${Ca}$ increases along the black arrow, starting at 
${Ca} = 0.05$ with increments of 
$0.05$. Windows of 
${\pm }5$ nm around the steady mean value are shown for (b) 
${Ca} = 0.20$ and (c) 
${Ca} = 0.25$. With red arrows we indicate thermal oscillations and motion similar to stick-slip.
Overview of MD simulations carried out in this work. For each simulation, 
${Ca}$ and steady 
$\Delta x$ (in nm) are given. If simulation is unstable, then ‘unst.’ is reported instead of a 
$\Delta x$ value.
$^a$Calibration simulation.
$^{b}$Stick-slip like behaviour.
The largest stable 
${Ca} = 0.25$ at 
$\theta _0=95^\circ$ exhibits oscillations similar to the so-called ‘stick-slip’ behaviour (Orejon, Sefiane & Shanahan Reference Orejon, Sefiane and Shanahan2011; Varma, Roy & Puthenveettil Reference Varma, Roy and Puthenveettil2021). 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 
$\Delta 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}, \theta _0) = (0.25,90^\circ )$ is distinct from the oscillations observed for 
${Ca} = 0.20$. We show the enlarged view of 
$\Delta 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 the 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 
$\theta _0 = 127^\circ$ at 
${Ca} = 0.90$ and 
$\theta _0 = 38^\circ$ at 
${Ca} = 0.02$. Simulations exhibiting stick-slip oscillations are indicated with a superscript 
$b$ in table 1.
3.3. Droplet break-up
As 
${Ca}$ is increased further, 
$\Delta x(t)$ measured from MD does not stabilize around some finite value. Instead, 
$\Delta x(t)$ grows continuously. An example of 
$\Delta x(t)$ at 
$({Ca},\theta _0) = (0.30,95^\circ )$ 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 
$\Delta 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 
$2 U_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 Reference Eggers2005) or hydrodynamic assist (Afkhami et al. Reference Afkhami, Buongiorno, Guion, Popinet, Saade, Scardovelli and Zaleski2018; Liu, Carvalho & Kumar Reference Liu, Carvalho and Kumar2019a; Fullana, Zaleski & Popinet Reference Fullana, Zaleski and Popinet2020).
Drop displacement in MD with 
$\theta _0 = 95^\circ$ at 
${Ca} = 0.30$. In insets (a–c), we show the drop shape at three selected time instances. The selected time is shown by a green circle with an arrow pointing to the inset. The green dashed line corresponds to the relative wall speed 
$2 U_w$.
4. Calibrating CFM models against MD
The aim of this section is to identify the free parameters in PF and VOF such that the displacement 
$\Delta x$ obtained from the continuum models matches the displacement obtained from MD simulations. The four calibration pairs 
$({Ca},\theta _0)$ are (marked with superscript 
$a$ in table 1) 
$(0.60,127^\circ )$, 
$(0.20,95^\circ )$, 
$(0.10,69^\circ )$ and 
$(0.015,38^\circ )$. For each 
$\theta _0$, we have chosen the largest steady 
${Ca}$ value available from MD simulations. If the chosen steady 
${Ca}$ value leads to large stick-slip like oscillations of 
$\Delta x$ (such as observed in figure 6c), then we select the previous (smaller) steady 
${Ca}$ value. The unsteady simulations are not suitable for calibration, because 
$\Delta 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 
$\rho _y$ over the height of the channel:
\begin{equation} \rho_\ell = \int_{y_b}^{y_t} \rho_y (y)\, {{\rm d} y} \approx 990\ \mathrm{kg}\ \mathrm{m}^{{-}3}. \end{equation}
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. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020). 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 
$\rho _v^{VOF} = 9.9\ {\rm kg}\ {\rm 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 (the distance between periodic boundary conditions) is matched with the distance between the left and right periodic boundary conditions in the MD system. The channel height is based on the chosen hydrodynamic wall position at the centre of the red bin (figures 5a–d), which results in 
$H = 29.22$ nm. This position is the subject of an investigation in itself (Herrero et al. Reference Herrero, Omori, Yamaguchi and Joly2019) and is discussed further in Appendix E.
The slip length (2.4) used in CFM models in principle can be related directly to MD simulations. Previous work (Huang et al. Reference Huang, Sendner, Horinek, Netz and Bocquet2008) has demonstrated that the slip length of the MD system follows a quasi-universal relationship with respect to equilibrium contact angle 
$\theta _0$. However, for the chosen MD system (SPC/E water on the 
${\rm SiO}_2$ surrogate wall), accurate slip length quantification has not yet been done. To obtain an indication about the validity of 
$\ell _s = 0$ for the selected MD system, we compare streamwise velocities near the wall between PF and MD (Appendix E). Through this comparison, we found that 
$\ell _s = 0$ holds for 
$\theta _0 = 38^\circ \text {--} 95^\circ$, while for 
$\theta _0 = 127^\circ$, the appropriate choice is 
$\ell _s = 0.44$ nm. These values are used in both VOF and PF as input parameters without any additional fitting; see the sixth and ninth columns of table 2. Note that 
$\ell _s$ in VOF is fixed and distinct from 
$\lambda$ – the length scale in condition (2.8), which is used for calibration.
Results of CFM calibration against MD. For each 
$({Ca},\theta _0)$ pair, we report MD reference data (steady displacement 
$\Delta x$) in the third column. In the following four columns, we report VOF parameters (2.4) and (2.8), and resulting steady displacements. In the final five columns, PF parameters (2.4), (2.9) and (2.12), and resulting displacements, are given.
4.1. VOF calibration
Figure 8 compares 
$\Delta x$ obtained from VOF (red) and MD (black) for the four calibration configurations. We have adjusted the parameter 
$\lambda$ in (2.8) to find the best fit of 
$\Delta x$ as 
$t\rightarrow \infty$. We observe that the droplet displacement from VOF stabilizes at values that agree well with the atomistic simulations, with the most challenging calibration configuration being 
$\theta _0=127^\circ$ (figure 8a). The 
$\lambda$ values that reproduce the displacement obtained from MD are listed in the sixth column of table 2. For the hydrophobic configurations, we obtained 
$\lambda = 4.325$ nm for 
$\theta _0 = 127^\circ$ and 
$\lambda = 0.935$ nm for 
$\theta _0 = 95^\circ$. We observe that for 
$\theta _0 = 95^\circ$, 
$\lambda$ is roughly four times smaller than for 
$\theta _0 = 127^\circ$. Naively, 
$\lambda$ can be regarded as a slip-related length scale near the contact line. Larger 
$\lambda$ for 
$\theta _0 = 127^\circ$ then suggest smaller friction, in line with the findings from PF calibration (§ 4.2).
Time evolution of 
$\Delta x$ in MD, PF and VOF for all calibration configurations.
For 
$\theta _0 = 69^\circ$ and 
$38^\circ$, we found even smaller 
$\lambda$ values, 
$0.313$ nm and 
$0.146$ nm, respectively. Qualitatively, this continuous behaviour of decaying 
$\lambda$ for smaller 
$\theta _0$ follows the contact line friction argument presented later (§ 4.2). It is also interesting to note that the obtained 
$\lambda$ values from this calibration procedure are of order 0.1–4 nm. This is similar to what is used by Legendre & Maglio (Reference Legendre and Maglio2015), who set 
$\lambda = 1$ nm in macroscopic simulations. The mesh spacing in the VOF simulations is 
$\varDelta /2 = 0.457$ nm. For 
$\theta _0 = 38^\circ$ and 
$69^\circ$, we have 
$\lambda < \varDelta /2$ (table 2). On the other hand, for 
$\theta _0 = 95^\circ$ and 
$127^\circ$, we have 
$\lambda > \varDelta /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 
$\theta _{num}$ for 
$\theta _0 = 95^\circ$ and 
$127^\circ$ should be viewed as a numerical parameter to tweak the curvature near the wall.
4.2. PF calibration
The droplet displacement 
$\Delta 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.
4.2.1. PF calibration proposed in the literature
For the Cahn–Hilliard PF model, there is a standard calibration procedure proposed by Yue et al. (Reference Yue, Zhou and Feng2010) and Yue & Feng (Reference Yue and Feng2011). The sequential steps are: (i) choosing the interface thickness 
$\epsilon$ that is suitable to describe the physical problem; (ii) setting the PF mobility 
$M$ according to the sharp interface limit (Yue et al. Reference Yue, Zhou and Feng2010); and (iii) calibrating the contact line friction 
$\mu _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 \approx 38$ nm (figure 1a). Based on this, the interface thickness is set to 
$\epsilon = 0.7$ nm. According to the sharp interface limit (Yue et al. Reference Yue, Zhou and Feng2010), the criterion for choosing 
$M$ is
\begin{equation} M > 1/16 \epsilon^2 / \sqrt{\mu_v \mu_\ell} . \end{equation}
Inserting the chosen values for interface thickness, vapour viscosity and liquid viscosity, we obtain 
$M > 3.21 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1}\ {\rm s}^{-1}$. In the MD simulations, we do not observe any physical effect that would hint towards a large 
$M$ value. Therefore, for all 
$\theta _0$ in this section, we select 
$M = 3.5 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1} \ {\rm s}^{-1}$, which is close to the lower limit.
The last step is to carry out simulations with different 
$\mu _f$ until a steady state is attained that matches the displacement obtained from MD. Table 3 summarizes the displacement obtained for each calibration pair 
$({Ca},\theta _0)$ at 
$\mu _f = 0, 0.1, 0.2, 0.3, 2.36,11.84$. The fourth column of table 3 shows the displacements for 
$\mu _f = 0$. We observe that no steady solution has been obtained for 
$\theta _0 = 127^\circ$. The steady 
$\Delta x$ obtained for 
$\theta _0=95^\circ$ overestimates the MD value (third column of table 3) by roughly 
$50\,\%$. Increasing 
$\mu _f$ only deteriorates the agreement. Therefore, we conclude that the matching procedure proposed by Yue & Feng (Reference Yue and Feng2011) is not adequate to calibrate the PF model for the chosen nanoscale configuration at 
$\theta _0 = 127^\circ$ and 
$95^\circ$.
Calibrating PF with MD by changing 
$\mu _f$ (strategy proposed by Yue & Feng Reference Yue and Feng2011). In the third column, we show the reference 
$\Delta x$ from MD. To the right, we show steady 
$\Delta x$ (in nm) obtained from PF for specified 
$\mu _f / \mu _\ell$. If no steady state exists, then we write ‘unst.’. In all simulations, 
$M = 3.5 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1}\ {\rm s}^{-1}$ and 
$\epsilon = 0.7$ nm.
For 
$\theta _0 = 69^\circ$ and 
$38^\circ$, on the other hand, 
$\mu _f = 0$ results in an underestimated drop displacement (table 3) compared to MD. We observe that the required PF contact line friction for 
$\theta _0 = 38^\circ$ (
$\mu _f = 11.84\mu _\ell$) is roughly five times larger than the one for 
$\theta _0 = 69^\circ$ (
$\mu _f = 2.36\mu _\ell$). More hydrophilic 
$\theta _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 
$\mu _f$ values.
4.2.2. Calibrating PF by adjusting mobility for 
$\theta _0 = 127^\circ$ and $95^\circ$
A way to reduce friction near the contact line is to increase 
$M$ (figure 3a) and thus allow for more diffusion. Therefore, we impose 
$\mu _f = 0$ and increase 
$M$ for 
$\theta _0 = 127^\circ$ and 
$95^\circ$ until the steady 
$\Delta x$ from the PF agrees with MD. The required PF mobility values are reported in table 2. We observe that for 
$\theta _0 = 95^\circ$, a three times larger 
$M$ is required compared to the 
$M$ set by the sharp interface limit (see the tenth column of table 2), whereas for 
$\theta _0 = 127^\circ$, 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 
$\epsilon < 4\,M^{1/2}\mu _v^{1/4}\mu _\ell ^{1/4}$. As we increased 
$M$, the interface thickness 
$\epsilon = 0.7$ nm was kept constant, therefore the condition is satisfied. However, this calibration procedure produces a noticeable diffusion near the contact line (figure 23c), which is not observed in the MD simulations.
5. 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 behaviour of the droplet. Moreover, an important practical aspect is the accuracy of the CFM when it comes to predicting the droplet behaviour to characterize the sheared droplet system quantitatively for a range of capillary numbers. Specifically, in this section, we fix the parameters for PF and VOF to values reported in table 2.
5.1. Time evolution of drop displacement
Figure 9 shows the droplet displacement as a function of time for a fixed 
$\theta _0=95^\circ$ but different 
${Ca}$. Figure 9(d) is identical to figure 8(b) and corresponds to the conditions for which the system was calibrated, 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 
$\Delta x$ at the same time instant) than VOF and MD.
Time evolution of 
$\Delta x$ in MD, PF and VOF for 
${Ca} \in (0.05, 0.30)$, with equilibrium contact angle 
$\theta _0 = 95^\circ$. The calibrated functions 
$\Delta x (t)$ from PF and VOF are shown in (d) with green and red lines, respectively. The a priori measurements of 
$\Delta x (t)$ from MD (§ 3) are shown with a black line. This colour code is retained for all comparisons that follow.
Figures 9(a,b,c,e) show the transient behaviour of CFM models in off-calibration conditions for 
${Ca} = 0.05$, 
$0.10$, 
$0.15$ and 
$0.25$. For moderate capillary numbers 
${Ca} \geqslant 0.10$, we observe qualitatively similar behaviour as for 
${Ca} = 0.20$. The CFM models predict the transient and steady 
$\Delta x$ values rather accurately. The PF, however, is always slightly slower compared to the VOF model. For 
${Ca} = 0.05$, the steady 
$\Delta 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.
To conclude the 
$\theta _0 = 95^\circ$ 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 
$\Delta x$ observed from MD. This is in line with observations in steady situations (figures 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 
$\theta _0=38^\circ, 95^\circ$ and 
$127^\circ$. Qualitatively similar results to those observed in figure 9 are obtained; see supplementary figures 1–3 available at https://doi.org/10.1017/jfm.2022.219.
5.2. Interface shape
For interface shape comparisons, the data are presented as the variation of the angle along the interface 
$\theta ( y )$ (figure 2a). Figure 10 compares 
$\theta (y)$ obtained from MD, PF and VOF for 
$\theta _0=95^\circ$ and different 
${Ca}$. The equilibrium angles at the bottom and top walls are shown with black dotted lines, while the hydrodynamic wall positions are presented with black dashed lines. The interface shapes in both calibration (
${Ca}=0.20$) and in off-calibration conditions (
${Ca} = 0.05$, 
$0.10$, 
$0.15$ and 
$0.25$) are similar between the models. We observe that both PF and VOF exhibit more pronounced differences from MD at the top wall (near the advancing contact line). At 
${Ca} = 0.25$ (figure 10e), the CFM model predictions are notably different from MD data compared to the other capillary numbers. The loss of accuracy could arise from the large 
$\Delta x$ oscillations in MD at 
${Ca} = 0.25$ (§§ 3 and 6.2) that are not captured with the CFM models.
Steady interface shape from MD, PF, and VOF simulations, 
$\theta _0 = 95^\circ$. Equilibrium angles and solid wall locations are shown with black dotted and dashed lines, respectively.
The MD data in figure 10 do not extend all the way to the wall, so the exact value of the dynamic contact angle is not known. At 
$\theta _0=95^\circ$, the dynamic contact angle in PF simulations is equal to the equilibrium angle since 
$\mu _f = 0$ (table 2). For VOF, the dynamic contact angle is different from 
$\theta _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 
$\lambda = 0.935 \ {\rm nm} > 0.467\ {\rm nm} = \varDelta / 2$ (table 2).
We have repeated this comparison of interface shapes at different 
${Ca}$ values between MD, PF and VOF for equilibrium contact angles 
$\theta _0 = 127^\circ$, 
$69^\circ$ and 
$38^\circ$. The results are reported in supplementary figures 4–6 and are similar to what is presented above.
5.3. Steady drop displacement
Figure 11 shows the steady displacement as function of 
$Ca$ for different 
$\theta _0$. As expected, the MD, PF and VOF points collapse for the calibration configurations (marked with vertical arrows). For the 
$\theta _0 = 127^\circ$ configuration, the PF and VOF models overestimate the steady displacement for large 
${Ca}$. A possible reason behind the discrepancy between PF, VOF and MD could be the slip length 
$\ell _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 Reference Thompson and Troian1997). This, consequently, would reduce 
$\Delta x$ and possibly move the PF and VOF predictions closer to the MD results. For 
$\theta _0 = 38^\circ$, the agreement between PF, VOF and MD diverges as we increase 
${Ca} > 0.02$. A reason for the disagreement could be the increased contact line friction (§ 6.1) at receding contact line, from which a liquid film is formed (§ 5.5 and figure 14). This asymmetry is not taken into account in the CFM simulations.
Steady 
$\Delta x$ as a function of 
${Ca}$ for all 
$\theta _0$.
5.4. Critical capillary number
Figure 12 shows the critical capillary number, 
${Ca}_c$, as a function of 
$\theta _0$ from the three methods. To determine 
${Ca}_c$ for a given 
$\theta _0$, we take the mean of the largest steady (
${Ca}_s$) and smallest unsteady (
${Ca}_u$) capillary numbers that were simulated, i.e.
\begin{equation} {Ca}_c = \frac{{Ca}_u + {Ca}_s}{2} \pm \frac{{Ca}_u - {Ca}_s}{2}. \end{equation}
The uncertainty is determined as half of the difference between these capillary numbers. To reduce 
$\Delta {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 
$\theta _0$, which has also been reported in previous numerical (Sbragaglia et al. Reference Sbragaglia, Sugiyama and Biferale2008) and analytical (Hocking Reference Hocking2001; Eggers Reference Eggers2004) investigations.
Critical capillary number 
${Ca}_c$ as a function of 
$\theta _0$.
We did not manage to run any VOF simulations for 
$\theta _0 = 127^\circ$ and 
${Ca} \geqslant 0.80$. Due to a large 
$\lambda$ (table 2), the dynamic angle for 
${Ca} \geqslant 0.80$ falls outside the physically admissible range 
$(0^\circ, 180^\circ )$. Consequently, in figure 12 we compare only PF predictions with the MD data at 
$\theta _0 = 127^\circ$. The PF model slightly under-predicts the 
${Ca}_c$ value. This discrepancy could again be due to the uncertainty in the slip length 
$\ell _s$.
Figure 12 shows that predictions of 
${Ca}_c$ for 
$95^\circ$ and 
$69^\circ$ agree with MD results within the accuracy bounds of MD. For 
$\theta _0=38^\circ$, however, discrepancies are observed, where 
${Ca}_c$ values computed from both CFM models are larger compared to the MD results. Note that the differences at 
$38^\circ$ are enhanced due to the logarithmic scale of the 
${Ca}_c$ axis in figure 12.
5.5. Above the critical capillary number
In this subsection, 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, 
$\theta _0 = 95^{\circ }$ at 
${Ca} = 0.30$, and 
$\theta _0 = 38^{\circ }$ at 
${Ca} = 0.05$.
Figure 13 shows the drop shape at four time instances for 
$({Ca},\theta _0)=(0.30,95^{\circ })$. We observe that the three models provide similar deformed states at 
$t \approx 5$ ns and 
$t \approx 10$ ns; see the top two rows in figure 13. The time instance just before and right after the break-up is shown in the third and fourth rows of figure 13. The thread connecting the lower and upper drops is very thin in the MD and VOF simulations, compared to PF. 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 a slightly thinner neck (figure 13c) and VOF with a very small drop neck (figure 13k). The time instants at which the break-up occurs are remarkably similar: 
$t^{MD}_b = 23.41$ ns for MD, 
$t^{PF}_b = 24.37$ ns for PF, and 
$t^{VF}_b = 23.23$ ns for VOF. For a complete time-dependent animation of MD, PF and VOF simulations side by side, see supplementary movie 1.
Snapshots of interface shape evolution over time from MD (a–d), PF (e–h) and VOF (i–l), with equilibrium contact angle 
$\theta _0 = 95^\circ$ and capillary number 
${Ca} = 0.30$.
We turn our attention to 
$\theta _0 = 38^{\circ }$ and 
${Ca} = 0.05$. To avoid a premature break-up of the liquid bridge in the PF simulation, we use 
$\epsilon = 0.35$ nm for this particular simulation since in this unsteady example, a third length scale – width of the liquid bridge in the deformed state – becomes important. Droplet shapes at five time instances from MD, PF and VOF are shown in figure 14. We observe that initially (top two rows in figure 14) the drop is deforming only slightly. After a certain time (figures 14c,h,m), a liquid film is generated at the receding contact lines of the drop. Very similar drop shapes are observed in MD, PF and VOF at different time instances. The VOF simulation takes roughly 
$5$ ns less to form a film similar to one observed in MD. For PF, one has to wait for around 
$15$ ns more than in MD. As time progresses, the liquid film becomes longer (figures 14d,i,n), and the liquid bridge becomes thinner. We also note that the tip of the film forms a thicker drop-like region. Finally, the liquid film from the periodic image merges with the liquid bridge (figures 14e,j,o). For MD and VOF simulations, the coalescence occurs only near one of the walls. For PF, on the other hand, a symmetric configuration is obtained with fully wetted top and bottom walls. Animations of the simulations can be found in supplementary movie 2.
Snapshots of interface shape evolution over time from MD (a–e), PF ( f–j), and VOF (k–o), with equilibrium contact angle 
$\theta _0 = 38^\circ$ and capillary number 
${Ca} = 0.05$. For PF, we have 
$\epsilon = 0.35$ nm.
6. Molecular physics of the sheared droplet
In this section, we present molecular phenomena of the sheared droplet that are particularly challenging to deal with 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.
6.1. Contact line friction measurements from MD
To extract the 
$\mu _f$ from MD, we use the formula proposed by Yue & Feng (Reference Yue and Feng2011):
\begin{equation} \left[\frac{\sqrt{2}}{3}\,\frac{\mu_f}{\mu_\ell}\sin\theta\right]{Ca} = \cos\theta_0 - \cos\theta . \end{equation}
Here, 
$\theta$ 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, Blake & De Coninck (Reference Fernández-Toledano, Blake and De Coninck2019) and Fernández-Toledano et al. (Reference Fernández-Toledano, Blake, De Coninck and Kanduč2020). For this work, however, we have determined that a non-equilibrium approach based on fitting (6.1) to MD data is the most efficient approach.
To determine the 
$\mu _f$ from (6.1), the dynamic contact angle 
$\theta$ has to be determined first for each 
${Ca}$ value. We extract the dynamic contact angle above the first reliable bin (green bins in figures 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 
$\theta _0$ at the same distance from the wall. The obtained dynamic contact angles for all 
$\theta _0$ values are gathered in figure 15. For each 
$\theta$, we also display error bars, obtained by adding 
$\pm$2 to the polynomial order of the interpolation. As expected, when 
${Ca}$ increases, so does the deviation of 
$\theta$ from 
$\theta _0$.
Measured dynamic contact angles 
$\theta$ for different 
${Ca}$ along with uncertainty intervals for advancing (red squares) and receding (blue rhombus) contact lines. The dashed lines are the best least-squares fits of (6.1). The equilibrium contact angles are (a) 
$\theta _0 = 127^\circ$, (b) 
$95^\circ$, (c) 
$69^\circ$, and (d) 
$38^\circ$.
We use a least-squares fit to match (6.1) to 
$\theta$ measurements. The fit provides 
$\mu _f/\mu _\ell$ 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 frictions are listed in table 4 along with previously reported 
$\mu _f$ from calibration of PF simulations.
Comparison of contact line friction used in PF simulations and values obtained directly from MD results. For PF, the same contact line friction is used for advancing and receding contact lines, while in MD, the contact line friction can be different.
We observe that larger line friction parameters are measured for smaller 
$\theta _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 Reference Blake and Haynes1969). Interestingly, we observe a difference between advancing and receding line friction for 
$\theta _0 = 127^\circ$ and 
$38^\circ$. Future 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 linearized MKT model – where line friction is defined in a slightly different way – would not change the conclusions. The line frictions 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 
$\mu _f$ obtained from MD directly in the PF boundary condition (2.12).
6.2. Stick-slip like oscillations
To differentiate between fluctuations caused by molecular-scale motion and large stick-slip 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 
$\sigma l_{th}^2$, giving 
$l_{th}=\sqrt {k_B T/\sigma }\simeq 0.27$ nm. The other two scales are the van der Waals diameter of SPC/E water molecules 
$\sigma _{SPC/E} \simeq 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 
${\rm RMSF}=(\langle \Delta x^2\rangle -\langle {\Delta x}\rangle ^2)^{1/2}$ for each of the four contact lines, both at equilibrium and under shear conditions. For 
$\theta _0 = 95^\circ$, we show the obtained RMSF values as a function of 
${Ca}$ in figure 16(a).
(a) Root-mean-square fluctuations (RMSF) of all four contact lines against capillary number for 
$\theta _0=95^\circ$. The error bars show a conservative estimate of the standard error. The spacing in the hexagonal silica lattice (blue dashed line), the van der Waals (vdW) diameter of the SPC/E water molecules (red dot-dashed line), and the length scale of thermal fluctuations (green dotted line) are shown. (b) Drop displacement over time for 
$\theta _0=95^\circ$ and 
${Ca}=0.25$. The mean (
$\Delta x$), left (
$\Delta x_l$) and right (
$\Delta x_r$) displacements are shown. A close-up on a time interval showing pinning/depinning is shown in the inset. The green dotted lines show the slope corresponding to wall velocity: at the beginning of the simulation, both advancing and receding contact lines stick to the wall and thus move apart from each other with velocity 
$2U_w$; during stick-slip events, the receding contact line sticks to the wall while the advancing line maintains a steady motion, thus the displacement matches the velocity of a single wall 
$U_w$.
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 (de-pinning) 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 that 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 more closely the MD simulation with 
$\theta _0=95^\circ$ 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, 
$\Delta x_l$ and 
$\Delta x_r$ are synchronized in an oscillatory motion. However, there are a few time intervals (between 17 and 20 ns, and between 45 and 47 ns) where indeed complete stick-slip occurs. In these intervals, the advancing and receding speeds match the magnitude of wall velocity; see the inset of figure 16(b).
It appears that local (pinning/depinning) and global (oscillations) processes co-exist. 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 assumed implicitly 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 longer 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 Reference Hocking2001; Eggers Reference Eggers2005; Varma et al. Reference Varma, Roy and Puthenveettil2021). However, to the best of the authors’ knowledge, it has not been observed directly.
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 
$\theta _0$. The distribution of the contact angle oscillations has been identified previously (Smith et al. Reference Smith, Müller, Craster and Matar2016). Fluctuations of the contact angle would in turn induce oscillations in 
$\Delta x (t)$.
6.3. Three-dimensional nature of drop break-up
In figures 13(d,h), we have observed that the interface shapes obtained from MD and PF display more diffuse regions at the tip of the broken thread. It is tempting to conclude that PF captures the MD behaviour correctly. However, the MD snapshots in figure 13 have been averaged in the spanwise 
$z$ direction, while PF simulation is a purely two-dimensional (2-D) 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 
$\theta _0 = 95^\circ$ 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 (3-D) holes (figure 17b) before disconnecting completely. Averaging in the 
$z$ direction results in a lower density (than the bulk), thus giving the impression of a diffuse interface (figure 13e).
Detailed views of the molecular system upon breakage for 
$\theta _0=95^\circ$ and 
${Ca}=0.4$, observed at 
$t=9.45$ ns: (a) overview from the side, and (b) close-up on the neck region. Two transparent periodic images are added at the sides. The close-up (b) is obtained by positioning the camera orthogonal to the thin thread in (a).
It is also worth noticing that the formation of these 3-D threads occurs very quickly. The time that 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 2-D 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 2-D and 3-D breakage dynamics. This in turn suggests that the selected MD system dimensions are reasonable to produce 2-D results and can be compared directly with 2-D CFM simulations.
7. 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 modelling the molecular physics in a nanoscale channel.
7.1. Challenges of the chosen CFM models
For 
$\theta _0 = 127^\circ$, 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 that we have chosen for calibration. For VOF, it is important to consider condition (2.8) as a numerical model that allows adjustment of the interface curvature at the contact line. For 
$\theta _0 = 127^\circ$, this was achieved by increasing 
$\lambda$ to values larger than the grid size, i.e. 
$\lambda > \varDelta /2$ (table 2). This leads to a negative 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(\theta _{num}) > G(\theta _0)$ and thus 
$\theta _{num} > \theta _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 therefore does 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 
$\theta _0$. Here, mobility is considered as a numerical tuning parameter to match the droplet displacement and does not correspond to the actual molecular diffusion at the interface (see Appendix F and figures 23a–c).
Very low contact angles impose other challenges for CFM. For 
$\theta _0 = 38^\circ$, we observed that steady drop displacement diverged (figure 11d) between the three models. Indeed, the prediction of a critical capillary number showed deviations between the methods (figure 12). In addition, the time that 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 
$\lambda$ obtained through VOF calibration (§ 4.1) is reasonable (
$\lambda < \varDelta / 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 
$\theta _0 = 38^\circ$, we observed significantly larger 
$\mu _f$ for the receding contact line compared to that for the advancing contact line. In the CFM models, on the other hand, the contact line properties (
$\mu _f$ for PF and 
$\lambda$ 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 for the CFM model inaccuracy.
Finally, there are limitations of the chosen CFM at very low and high 
${Ca}$ values. At high 
${Ca}$ values (close to 
${Ca}_c$), we observed enhanced oscillations of 
$\Delta x (t)$ (figures 6(c) and 16, superscript 
$b$ in table 1). More detailed analyses of these oscillations were presented in § 6.2 for 
$\theta _0 = 95^\circ$. The CFM model does not include the intrinsic oscillations present in the molecular reality. Consequently, the agreement between CFM and MD in drop time evolution (figure 9(e) and supplementary figure 6) is degraded. This could also be a potential source of increased discrepancy between interface shapes from CFM and MD (figure 10e). Also, for low 
${Ca}$ values (for example, 
${Ca} = 0.05$ in figure 9a), we have observed a 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}$ values. Consequently, the signal-to-noise ratio in MD is much smaller for smaller 
${Ca}$ values, and the large oscillations can give the impression of larger inaccuracy in CFM models.
7.2. 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 
$\ell _s$ exhibits only small variations in 
$\theta _0$. Below bulk liquid, the slip length had to be 
$0.44$ nm for 
$\theta _0 = 127^\circ$, and 
$0$ nm for all other 
$\theta _0$. As discussed in § 7.1, for low-friction configurations, 
$\theta _0>90^\circ$, it was necessary to adjust 
$\lambda$ 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 
$\ell _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. 
$\theta _{rec}\sim 0$ for 
${Ca}\ll 1$). On the other hand, the interpretation for the almost vanishing receding friction on the 
$\theta _0 = 127^\circ$ surface is less obvious. To the best of the authors’ knowledge, nowhere before has an asymmetric behaviour 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 hydrophobic substrates are low-friction surfaces – under doubt. Indeed, it is not clear whether the value of contact line friction can be predicted from 
$\theta _0$ alone (Liu, Yu & Wang Reference Liu, Yu and Wang2019b; Wang Reference Wang2019). Reasoning with the frame of mind of KMT instead, line friction asymmetry hints at 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 
$\mu _f$ values in PF and different 
$\lambda$ values in VOF for advancing and receding contact lines.
7.3. Potential modelling directions
It has been proposed previously 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 (Reference Blake1993). This approach was later evaluated against MD simulations (Qian et al. Reference Qian, Wang and Sheng2003; Qian, Wang & Sheng Reference Qian, Wang and Sheng2006; Mohand et al. Reference Mohand, Hoang, Galliero and Legendre2019), and good agreement has been found. However, recently Lācis et al. (Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020) 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 (
$\ell _s = 0.44$ nm for 
$\theta _0 = 127^\circ$, and 
$\ell _s = 0.0$ nm for 
$\theta _0 = 38^\circ \text {--}95^\circ$). Previously (Lācis et al. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020), for 
$\theta _0 = 95^\circ$ we extracted 
$\ell _s = 0.17$ nm, which is close to what we have in the current work. However, the most important overlap between our current and previous (Lācis et al. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020) 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. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020) 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. Reference Hartmann, Fricke, Weimar, Gründing, Marić, Bothe and Hardt2021). Another effect, which was not considered in this work, is so-called disjoining pressure (Pismen & Pomeau Reference Pismen and Pomeau2000). Including it would also allow for direct modification of contact line motion. Furthermore, more detailed studies of local rheological effects (including the 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. An infinite friction parameter renders GNBC ineffective, whereas gradually reducing the 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 Reference Hadjiconstantinou1999; Zhang, Borg & Reese Reference Zhang, Borg and Reese2017; Borg et al. Reference Borg, Lockerby, Ritos and Reese2018) 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 Waals (Laurila et al. Reference Laurila, Carlson, Do-Quang, Ala-Nissila and Amberg2012) or Cahn–Allen (Eggleston, McFadden & Voorhees Reference Eggleston, McFadden and Voorhees2001). The level-set model (Tornberg & Engquist Reference Tornberg and Engquist2000) or lattice Boltzman method (Chen et al. Reference Chen, Kang, Mu, He and Tao2014) are other potential candidates for the simulation of multiphase flows; for instance, Latva-Kokko & Rothman (Reference Latva-Kokko and Rothman2007) showed how a no-slip LBM model is able to capture automatically 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, Roccon & Soldati (Reference Soligo, Roccon and Soldati2021). It must be recognized that the hybrid models (Zhang et al. Reference Zhang, Borg and Reese2017; Liu et al. Reference Liu, Zhang, Capobianchi, Borg, Zhang and Wen2021) 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 stick-slip like oscillations of the contact line are important, then some other means of continuum modelling can be pursued. For example, the fluctuating hydrodynamic interfaces model proposed by Flekkøy & Rothman (Reference Flekkøy and Rothman1996) and Smith et al. (Reference Smith, Müller, Craster and Matar2016) could be evaluated as a suitable choice. Finally, we mention the possibility of regularizing the contact line stress singularity via Brinkman's model for porous surfaces (Devauchelle, Josserand & Zaleski Reference Devauchelle, Josserand and Zaleski2007).
Despite the missing physical mechanisms in PF and 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 
$\theta _0$.
8. Conclusions
We have calibrated a standard Cahn–Hilliard PF model, as well as a standard geometric volume-of-fluid model with Cox-like wetting condition, against MD simulations of water over a no-slip substrate. The no-slip behaviour in the MD system is an outcome of electrostatic bonds between polar water molecules and polar wall molecules. Two parameters (mobility 
$M$ and contact line friction 
$\mu _f$) were adjusted in PF simulations, and one parameter (Cox cut-off length scale 
$\lambda$) in VOF simulations. Four different equilibrium contact angles (
$\theta _0 = 127^\circ$, 
$95^\circ$, 
$69^\circ$ and 
$38^\circ$) were investigated. For each 
$\theta _0$, the largest steady stick-slip free simulation was selected for calibration. The PF and VOF models were calibrated to match the steady 
$\Delta x$ – single scalar macroscopic measurement – observed in MD. Using the calibrated parameters, a series of simulations was 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 
$\theta _0 = 95^\circ$ and 
$\theta _0 = 38^\circ$. The CFM models predicted all qualitative features of the MD simulations. For 
$\theta _0 = 95^\circ$, drop displacement and break-up in half were predicted. For 
$\theta _0 = 38^\circ$, 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 
$\mu _f$ values were asymmetric between advancing and receding contact lines for 
$\theta _0 = 127^\circ$ and 
$38^\circ$.
Finally, we have discussed the variations of PF and VOF parameters for matching the MD results for all 
$\theta _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 
$\theta _0 \geqslant 95^\circ$. 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 receding contact lines is the source of the inaccuracies in PF and VOF when predicting the 
$\theta _0 = 38^\circ$ results.
We have demonstrated that by calibrating the CFM once for each 
$\theta _0$ by targeting a single global measure 
$\Delta x$, it is possible to obtain many accurate predictions of interface shape 
$\Delta 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 (
$\theta _0 = 95^\circ$ and 
$69^\circ$). On the other hand, more deviations were observed for hydrophobic (
$\theta _0 = 127^\circ$) and hydrophilic (
$\theta _0 = 38^\circ$) configurations. These accuracy limits must 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 material and movies
Supplementary figures, movies and data files for easier figure reproductions are available at https://doi.org/10.1017/jfm.2022.219.
Acknowledgements
Numerical simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC, NSC and HPC2N. The authors thank T. Fullana and Dr P. Johansson for fruitful discussions. The developers of the PARIS code are acknowledged for continuous advancement and support of the tool. S.Z. thanks the Dean of the Faculty of Sciences of Sorbonne Université, Stéphane Régnier, for the grant of a sabbatical leave that made his contribution to this work possible.
Funding
U.L., S.B., J.S., B.H. and M.P. were supported by Swedish Research Council (INTERFACE centre and grant no. VR-2014-5680). U.L. was supported by the ERDF (project no. 1.1.1.1/20/A/070). S.Z. acknowledges funding from ERC (advanced grant TRUFLOW).
Data availability statement
The dataset containing all the flow field maps used for the calibration of CFM methods is openly available under Creative Commons Attribution 4.0 International licence on Zenodo (Pellegrino & Hess Reference Pellegrino and Hess2022), at http://doi.org/10.5281/zenodo.5997091. Data files and scripts to plot figures in the main paper are provided in the supplementary material.
Declaration of interests
The authors report no conflict of interest.
Author contributions
U.L. and M.P. contributed equally to this work. G.A., S.B., B.H. and S.Z. conceived the original idea. U.L. performed the PF simulations, carried out MD post-processing and constructed PF/VOF/MD comparisons. M.P. configured the MD system, carried out MD simulations, improved flow field extraction from MD data, and extracted contact line friction from MD. J.S. performed the VOF simulations and implemented required customizations with feedback from S.Z. All authors analysed data. U.L., M.P., J.S. and S.B. wrote the paper with feedback from all co-authors.
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, customization for the 
$\theta _0 = 127^\circ$ configuration is explained.
A.1. Numerical implementation
All VOF simulations used a resolution of 
$N_x = 256$ (
$\varDelta _x = 0.624$ nm) cells in the streamwise direction, and 
$N_y = 32$ (
$\varDelta = \varDelta _y = 0.913$ nm) in the wall-normal direction. PARIS solves the general 3-D equations. The 2-D 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 built-in implementation of the algorithm by Weymouth & Yue (Reference Weymouth and Yue2010), 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, Zaleski & Bagheri (Reference Sundin, Zaleski and Bagheri2021). 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 Reference Afkhami and Bussmann2008). 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 
$\theta _0 = 95^\circ$, with 
${Ca} = 0.05$ and 
$0.15$. The time series of the drop displacements are presented in figure 18(a). For 
${Ca} = 0.05$, the drop displacement changed by 
$1.3\,\%$, and for 
${Ca} = 0.15$ it changed 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 Reference Legendre and Maglio2015). Nevertheless, we consider the observed convergence sufficient.
(a) VOF refinement study, showing drop displacement for 
$N_x = 256$ and 
$N_y = 32$ cells (solid lines) and 
$N_x = 512$ and 
$N_y = 64$ cells (dashed lines) for 
${Ca} = 0.05$ and 
$0.15$. (b) Curvatures estimated by the analytical relations (A2) and (A3) (solid and dashed lines, respectively).
A.2. Customisation for $\theta _0 = 127^\circ$
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 (
$|\boldsymbol {{f}}_\sigma |$, (2.6)). The expression for the curvature is
\begin{equation} \kappa = \frac{{\rm d}^2 h/{{\rm d} y}^2}{(1 + ({\rm d} h/ {{\rm d} y})^2)^{3/2}}, \end{equation}
where 
$h = h(y)$ is the height function giving the 
$x$ location of the interface. The sign of 
$\kappa$ also depends on the interface orientation; for brevity, this discussion is not included (see Aniszewski et al. Reference Aniszewski2021). 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 layers above the wall, respectively. From the definition of the angle, 
$h_0 = h_1 + \varDelta _y/\tan (\theta )$ (the sign in front of 
$\tan (\theta )$ depends on the interface orientation). We assume that 
$h_2 \approx h_1$. The curvature in the first cell layer becomes
\begin{equation} \kappa_{num} = \frac{(h_2 - 2h_1 + h_0)/\varDelta_y^2}{(1 + (h_2 - h_0)^2/(2\varDelta_y)^2)^{3/2} } \approx \frac{1}{\varDelta_y}\,\frac{1/\tan(\theta)}{(1 + 1/(4\tan^2(\theta)))^{3/2}}, \end{equation}
shown in figure 18(b) (solid line). The approximation of 
$\kappa _{num}$ is zero for 
$\theta = 0^\circ$, 
$90^\circ$ and 
$180\,^\circ$. As shown in the figure, 
$|\kappa _{num}|$ has one maximum in each of the intervals 
$0^\circ < \theta < 90^\circ$ and 
$90^\circ < \theta < 180^\circ$. For 
$\theta _0 = 127^\circ$, a minimum separation velocity was achieved for a steady-state receding contact line angle smaller than 
$180\,^\circ$, as expected from (A2). 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
\begin{align} \kappa_{num} = \frac{(h_2 - 2h_1 + h_0)/\varDelta_y^2}{(1 + ({-}h_3 + 6h_2 - 3h_1 - 2h_0)^2/(6\varDelta_y)^2)^{3/2} } \approx \frac{1}{\varDelta_y}\,\frac{1/\tan(\theta)}{(1 + 1/(9\tan^2(\theta)))^{3/2}}, \end{align}
where we assumed 
$h_3 \approx h_2 \approx h_1$. This expression results in significantly higher curvatures for extreme angles (figure 18b, 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 
$\theta _0 = 127^\circ$, not enough wall-normal heights could be computed at receding or advancing contact lines for curvature calculations. The current solution is therefore more robust. Wall-normal heights could be a viable solution if resolution is increased significantly.
Appendix B. Details of the Cahn–Hilliard PF model
In this appendix, we provide additional details of the Cahn–Hilliard PF model used in this work. The model is introduced briefly in § 2.2. The standard double-well potential is
\begin{equation} \varPsi (C) = \left( C+1 \right)^2\left( C-1 \right)^2/4. \end{equation}The wetting boundary condition (2.12) contains the so-called switch function, defined as
\begin{equation} g(C) = 0.5 - 0.75 C + 0.25 C^3. \end{equation}
This expression serves as a window function that ensures that the wetting boundary condition is applied only at the contact line. Furthermore, the cubic expression (B2) is not empirical. Instead, it is derived as the equilibrium solution of PF equations based on the selected potential (B1) 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 PF variable 
$C$ as
\begin{equation} \rho(C) = \rho_\ell\,\frac{C+1}{2} - \rho_v\, \frac{C-1}{2} \quad \mbox{and} \quad \mu(C) = \mu_\ell\,\frac{C+1}{2} - \mu_v\, \frac{C-1}{2}. \end{equation}
Recall that 
$\rho _\ell$ and 
$\mu _\ell$ are the density and viscosity of the liquid component, and 
$\rho _v$ and 
$\mu _v$ are the density and viscosity of the vapour component.
The introduced governing equations are linearized, written into the weak form, and solved using an open-source finite-element software FreeFEM (Hecht Reference Hecht2012), which allows easy specification of finite-element weak form. Linear elements were used for the PF 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 
$\Delta s_1 = 3.65$ nm far from the interface, and down to 
$\Delta s_2 = 0.24$ nm within the interface region. The constant time step was used through the simulation as 
$\Delta t = 0.002$ dimensionless time units. For very small 
${Ca}$, 
$\Delta t$ was reduced to ensure numerical stability. For simulations with smaller 
$\epsilon$, 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 the Github repository (Lācis & Bagheri Reference Lācis and Bagheri2020–2022). The data to reproduce the figures in the main paper are available in supplementary file 1.
Appendix C. Details of MD simulations
In this appendix, we provide the details necessary for the reader to understand the simulation procedure, interface extraction, and determination of equilibrium angle.
C.1. Numerical implementation
The thin quasi-2-D liquid meniscus (figure 1) is composed of 172 933 water molecules. The parameters for the SPC/E model are taken from the OPLS-AA force field. The parametrization of 
${\rm SiO}_2$ quadrupoles is summarized in table 5. Silicon atoms are treated as virtual sites without mass. Oxygen atoms are restrained to absolute coordinates by a spring of constant 
$\kappa _O$. The usage of position restraints grants the substrate some flexibility to reconfigure 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, 
$\delta 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. Reference Abraham, Murtola, Schulz, Páll, Smith, Hess and Lindahl2015).
Parametrization of the force field of silica quadrupoles (electrostatics excluded).
To obtain 
$\theta _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 Lennard-Jones (L-J) walls are placed beyond the silica surfaces at the location of periodic boundary condition. The L-J 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 19b). 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 wall and by 
$-d_{ne}$ on the bottom wall. The effective wall velocity is quantified as
\begin{equation} U_w = \frac{\delta x}{\delta t} = d_{ne}\,\frac{\delta\lambda}{\delta t} . \end{equation}
Here, 
$\delta \lambda$ is the increment of an auxiliary variable 
$\lambda \in [0,\infty )$ applied at each time step. The value 
$\lambda =0$ corresponds to configuration A, and 
$\lambda =1$ to configuration B. The desired wall velocity is obtained by setting the interpolation increment to 
$\delta \lambda =U_w\,\delta t/d_{ne}$.
Position restraints and wall treatment. (a) References for restrained positions are visualized via ‘ball-and-stick’ representation, while actual atoms of 
${\rm SiO}_2$ molecules are shown as transparent vdW spheres. (b) Auxiliary interpolation procedure for non-equilibrium simulations to produce an effective wall velocity 
$U_w$.
All simulations are performed in the NVT ensemble at 
$T=300\ {\rm K}$ and fixing the extent of the 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 L-J walls at 
$y=0$ and 
$y=L_y$ (figure 19a). A 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 lead to either local flow hindering or 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
\begin{equation} E_{kpp} = \tfrac{1}{2}m\left(\langle{c_x^2}\rangle+\langle c_y^2\rangle+\langle c_z^2\rangle\right) = \tfrac{3}{2}k_B T , \end{equation}
where 
$\boldsymbol {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
\begin{equation} \tfrac{3}{2}k_B T_0 = \tfrac{1}{2}m\left(u^2+\langle{c_x^2}\rangle+\langle{c_y^2}\rangle+\langle{c_z^2}\rangle\right). \end{equation}The difference between imposed and effective temperature is
\begin{equation} T_0-T = \frac{mu^2}{3k_B} \simeq \frac{m(\gamma/\mu)^2}{12k_B}\,{Ca}^2 = \varTheta\,{Ca}^2 > 0 , \end{equation}
where 
$\varTheta$ 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 
$\varTheta \simeq 0.78$ K, which entails that to cool down the near-wall molecules by 
$1$ K, one needs to prescribe at least 
${Ca}\simeq 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 coupling only the solid substrate and letting the liquid thermalize due to heat exchange. Trying this approach, we noticed that in the configurations with 
$\theta _0 \geqslant 95^\circ$ (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, Todd & Searles Reference Bernardi, Todd and Searles2010) or a dissipative-particle-dynamics thermostat (Soddemann, Dünweg & Kremer Reference Soddemann, Dünweg and Kremer2003; Goga et al. Reference Goga, Rzepiela, de Vries, Marrink and Berendsen2012). However, these thermostats are not currently implemented in GROMACS.
Hydrodynamic fields (density, velocity and temperature) are measured directly from MD trajectories. Each quantity is averaged in space on a grid with spacing 
$(h_x, h_y) \simeq (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 
$\rho ^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 initialization 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 potentially can hinder relaxation in equilibrium simulations and create velocity measurement artefacts in non-equilibrium ones. This averaging procedure is employed before interface extraction for both equilibration and sheared MD runs.
C.3. Interface extraction and $\theta _0$ measurement
In this subsection, we describe the extraction of the interface shape from the water density distribution 
$\rho (x,y)$. The distribution for the equilibration run with 
$\theta _0 = 95^\circ$ along the bin with vertical coordinate 
$y \approx 1.5$ nm is shown in figures 5(e, f). We consider two criteria to define the exact 
$x_{ni}$ coordinate of the interface, i.e.
\begin{equation} \rho(x_{ni}, y_n) = 0.5\rho_\ell \quad \mbox{and} \quad \rho(x_{ni}, y_n) = 0.5\rho_y (y_n). \end{equation}
These are based on global liquid density (C5a) and slice liquid density (C5b) at an 
$n$th vertical bin with coordinate 
$y_n$. If the sought density value is not located in any single bin, then linear interpolation is used to find the exact coordinate. The interface extraction according to (C5a) we call ‘global interface extraction’, and the obtained interface we call the ‘global interface’. This is the approach used typically in literature. However, the density layering (figure 20b) exhibits itself in the interface shape near the surface (figure 20a, green line). The interface extracted according to (C5b) we call the ‘slice interface’. The 
$x_{ni}$ obtained from (C5b) show reduced layering influence (figure 20a, blue line). Further away from the wall, results from (C5a,b) agree. From (C5b) it is also possible to define the interface point for adsorbed water layer (figure 20a, 
$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.
Equilibration MD run yielding 
$\theta _0 = 95^\circ$. Close-up near the bottom wall of extracted interface shape (a) and liquid density variation (b). Interface angle variation (c) along the curvilinear coordinate 
$s$ excluding the interface measurements closest to the walls. Inset (d) shows the interface angle along the vertical coordinate 
$y$ with all interface points included. Dashed lines in (c,d) correspond to the measured equilibrium contact angle 
$\theta _0$.
As the next step, we measure the equilibrium contact angle 
$\theta _0$. We use both interface shapes extracted according to (C5). From the interface shape, we compute the interface angle 
$\theta (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 in green in figure 20a). Therefore, angle measurements are located at the boundaries of the MD bins. The Young–Laplace equation for constant surface tension is
\begin{equation} \Delta p ={-} \sigma \gamma, \end{equation}
where 
$\Delta p$ is the pressure jump across the interface, and 
$\gamma$ 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 equation (C6), this results in constant curvature along the interface. We can express the interface angle as
\begin{equation} \theta(s) = c_1s + c_2, \end{equation}
where 
$c_1$ and 
$c_2$ are constants determined by boundary conditions, and 
$s$ is the curvilinear coordinate along the interface. For convenient comparison, we transfer 
$\theta (y)$ to 
$\theta (s)$. Since the theoretical function is known, we fit the obtained MD results with (C7). We observe that both global and slice interfaces deviate rapidly from linear relationships near the wall (figure 20d). The same effect is observed in figure 20(a), where the interface segments closest to the wall exhibit significantly different angles compared to segments 
${\approx }1.5$ nm above the wall. Therefore, we conclude that several interface points closest to the wall cannot 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 (C7) reaches its minimum. This process is applied to both global and slice interfaces. We 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 
$\theta (s)$ from the slice interface for the remaining interface points is presented in figure 20(c) with a blue line. The corresponding best linear fit is given with a red dotted line. To obtain the equilibrium angle 
$\theta _0$, we use the hydrodynamic wall position assumed in the main paper. The wall position is shown in red in figure 20(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 20c, red crosses). However, the total arc length of the interface depends on the contact angle. Therefore, the equilibrium angle 
$\theta _0$ is obtained iteratively in the following procedure. We first centre the MD data to ideal arc length for given 
$\theta _0$. Then the agreement with the given 
$\theta _0$ is verified by extrapolating the best linear fit of the MD data to wall locations. Finally, the next 
$\theta _0$ is set as the average of the previous estimate and the current estimate. The final obtained equilibrium angle in figure 20(c) is presented with black dotted lines. The inset (figure 20d) 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 cannot 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 
$\theta _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 (figures 5a–d), and the measured contact angles 
$\theta _0$ are 
$127^\circ$, 
$95^\circ$, 
$69^\circ$ and 
$38^\circ$, as stated in the main paper.
C.4. Polynomial extrapolation of MD interface shape
In this subsection, 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 
$\theta _0$ values (Appendix C.3), the first reliable bins for comparison with CFM have been identified (shown with green in figures 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 
$\theta _0 = 95^\circ$ in figure 21(a). With red dashed lines, we show the assumed wall positions. 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.
MD interface shape (a) and interface angle (b) for 
$\theta _0 = 95^\circ$ system at equilibrium and sheared (
${Ca} = 0.20$) configurations. (c) Convergence of 
$\Delta x$ with respect to polynomial fit order 
$N_p$. Wall locations are shown with red dashed lines in (a,b). The polynomial fit used to obtain the drop displacement estimate is shown with a purple dotted line (a). The equilibrium angle is shown with a black dotted line in (b). The purple crosses in (b) represent possible contact angle measurements for a few 
$N_p$ values.
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 21(c) for 
$\theta _0 = 95^\circ$ and 
${Ca} = 0.20$. Typically, we observe that initially there are large changes of the extrapolated drop displacement. However, after some order, the magnitude of the 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 
$\Delta 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 21b) and increase the order if the agreement is not satisfactory. The final chosen polynomial order for 
${Ca} = 0.20$ is 
$N_p = 7$. The final obtained steady displacement from MD is 
$\Delta x = 13.89$ nm. The shape of the fitted polynomial for the interface is shown in figure 21(a,b) with a purple dotted line.
While 
$\Delta 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 21(b) we add two more polynomial fits with 
$N_p = 8$ and 
$N_p = 11$. All three polynomial orders give very similar 
$\Delta x$ values (figure 21c). However, the predicted advancing dynamic contact angles at the wall (purple crosses in figure 21b) 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 21b). The polynomial extrapolation in the main paper is used to get a more accurate a posteriori steady 
$\Delta x$ measurement and to read off the dynamic contact angle at the reliable bin location. The convergence of 
$\Delta 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 gases) at 
$T = 300\ {\rm K}$ is
\begin{equation} p_{wv} = \frac{\exp(77.345 + 0.0057\,T - 7235/T)}{T^{8.2}} = 3523.88\ \mathrm{Pa}. \end{equation}The water vapour density is
\begin{equation} \rho_{wv} = 0.0022\,\frac{p_{wv}}{T} = 0.0258\ \mathrm{kg}\ \mathrm{m}^{{-}3}. \end{equation}
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\ {\rm K}$ gives us
\begin{equation} \mu_{wv} = \mu_{20^\circ} + \left(\mu_{50^\circ} - \mu_{20^\circ}\right) \frac{323.15\ {\rm K} - T}{30\ {\rm K}} = 1.04 \times 10^{{-}5}\ \mathrm{Pa} \ \mathrm{s}. \end{equation}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 show 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 22a). If the chosen wall location is appropriate, then 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 22a) 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.
(a) Assumed hydrodynamic wall position for the CFM simulations with respect to the molecular picture. (b–e) Horizontal velocity 
$u_x$ over horizontal slice at the bin roughly corresponding to the liquid density peak coloured green in (a). Results for 
$\theta _0 = 127^\circ$ (b), 
$95^\circ$ (c), 
$69^\circ$ (d) and 
$38^\circ$ (e) are reported. Calibration 
${Ca}$ numbers are stated in panel titles. Expected sensitivity to a slight shift in wall location is given as a grey shaded region.
We extract 
$u_x$ from MD at the green bin. The MD velocity data are obtained with the help of two averaging approaches. The first approach is the MD frame average over the steady regime together with a 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 right sides 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-averaged data until 
${\approx }10$ nm away from the contact line. The noisy data on the vapour side are 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 streamwise velocity distribution is shown in figures 22(b–e) with solid black lines. Contact angles and 
${Ca}$ numbers are presented in the titles of individual panels.
For comparison, we extract the 
$u_x$ distribution along the 
$x$ coordinate at the same 
$y$ location from all calibration no-slip PF simulations. The PF results for 
$\theta _0 = 127^\circ$ are shown in figure 22(b) 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 
$\ell _s$ value until a good match is obtained. Through this, we obtain 
$\ell _s = 0.44$ nm. The PF velocity results with 
$\ell _s = 0.44$ nm are shown in figure 22(b) with a solid red line. Now, a good agreement between MD and PF is obtained.
For 
$\theta _0 \leqslant 95^\circ$ (figures 22c–e), we observe that agreement between PF velocity predictions and MD results is good at 
$\ell _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 (
$\ell _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 areas (figures 22b–e) we show the regions in which the PF results would fall if the wall location were moved up or down by 
$0.1$ nm. These results are obtained by sampling 0.1 nm closer to 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 regions, we observe that the variations are very narrow for all 
$\theta _0$ values in the centre of the drop. However, approaching the contact line region for 
$\theta _0 = 95^\circ$ and 
$69^\circ$, 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 must be recognized that the slip length of MD systems has been studied extensively before; see, for example, works by Thompson & Troian (Reference Thompson and Troian1997) and Huang et al. (Reference Huang, Sendner, Horinek, Netz and Bocquet2008). In particular, Huang et al. (Reference Huang, Sendner, Horinek, Netz and Bocquet2008) investigated the slip length of SPC/E water over a range of surfaces, from silane monolayers to more common L-J models. They found that, up to a good accuracy, 
$\ell _s \sim ( 1 + \cos \theta _0 )^{-2}$. It must be noted that, to the best of authors’ knowledge, a similar study in MD with surfaces that form hydrogen bonds with water has not been carried out before, and it is out of scope of the current study. Nevertheless, the slip length values obtained in the current work agree qualitatively with results of Huang et al. (Reference Huang, Sendner, Horinek, Netz and Bocquet2008) – as 
$\theta _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 
$\theta _0 = 38^\circ \text {--}95^\circ$.
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 figures 5a–d). In addition, through determining the validity of the no-slip condition, we have demonstrated that PF can predict accurately 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 understanding of the influence of parameters 
$M$, 
$\mu _f$ and 
$\epsilon$ on the PF results, the steady flow field in the vicinity of the receding contact line is investigated. For this, the PF simulation with 
$\theta _0 = 95^\circ$ and 
${Ca} = 0.20$ is run until the steady 
$\Delta x$ is reached. The flow field at the last time instant is used to compute the streamlines. With black lines in figure 23, we present streamlines near the bottom left receding contact line in a zoomed-in window of roughly 
$12\ {\rm nm} \times 7\ {\rm 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.
Streamlines from PF simulations near the bottom left receding contact line. All simulations use equilibrium contact angle 
$\theta _0 = 95^\circ$ and capillary number 
${Ca} = 0.20$. The blue streamline describes a fluid parcel originating from within the liquid drop at 
$0.5$ nm distance from the wall (a second blue streamline has the same stream function value). The red isoline corresponds to the two-phase interface, defined as 
$C = 0$. With the light blue colour we show the variations of the 
$C$ function. PF mobility (a–c), contact line friction (d–f), and interface thickness (g–i) are varied. Other parameters in the corresponding row are kept constant. In (a–c), 
$\mu _f = 0$ and 
$\epsilon = 0.7$ nm. Then 
$\epsilon = 0.7$ nm and 
$M = 1.75 \times 10^{-15}\ {\rm m}^4\ {\rm N}^{-1} {\rm s}^{-1}$ are used in (d–f). In (g–i), we have 
$\mu _f = 0$ and 
$M = 1.08 \times 10^{-15}\ {\rm m}^4\ {\rm N}^{-1}\ {\rm s}^{-1}$.
In figures 23(a–c), the influence of PF mobility is shown. First, we investigate the flow field with the smallest 
$M$ (figure 23a). Overall, the streamlines are similar to the wedge solution derived and presented by Huh & Scriven (Reference Huh and Scriven1971), and the PF solution analysed thoroughly by Seppecher (Reference Seppecher1996). 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 23b), 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 a 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$ a hundred times larger (figure 23c), 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 recognized.
Figures 23(d–f) illustrate the effect of varying contact line friction 
$\mu _f$. When there is no friction (figure 23d), the contact angle at the wall is equal to 
$\theta _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 
$\mu _\ell$ (figure 23e) 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 23(d). The contact angle at the wall deviates from 
$\theta _0$. For the largest contact line friction (figure 23f), we observe an even more pronounced contact angle departure from 
$\theta _0$. In addition, the marked streamline crosses the interface and continues in the vapour phase.
Finally, in figures 23(g–i), we exemplify the effect of the interface thickness. As 
$\epsilon$ is reduced, the behaviour of the blue streamline changes from full crossover to vapour phase (figure 23i), to a small overshoot of an order of nm (figure 23h) and finally to no crossing over the interface (figure 23g). The interface shape, on the other hand, remains practically the same for 
$\epsilon = 0.70$ nm and 
$0.18$ nm, which is again a signature of the sharp interface limit.
Open access
