Steady moving contact line of water over a no-slip substrate

The movement of the triple contact line plays a crucial role in many applications such as ink-jet printing, liquid coating and drainage (imbibition) in porous media. To design accurate computational tools for these applications, predictive models of the moving contact line are needed. However, the basic mechanisms responsible for movement of the triple contact line are not well understood but still debated. We investigate the movement of the contact line between water, vapour and a silica-like solid surface under steady conditions in low capillary number regime. We use molecular dynamics (MD) with an atomistic water model to simulate a nanoscopic drop between two moving plates. We include hydrogen bonding between the water molecules and the solid substrate, which leads to a sub-molecular slip length. We benchmark two continuum methods, the Cahn--Hilliard phase field (PF) model and a volume of fluid (VOF) model, against MD results. We show that both continuum models can reproduce the statistical measures obtained from MD reasonably well, with a trade-off in accuracy. We demonstrate the importance of the phase field mobility parameter and the local slip length.


Introduction
The motion of a two-fluid interface contacting a flat solid surface poses a particularly difficult problem of continuum fluid mechanics. If the traditional point of view of a no-slip wall-a sharp transition between the phases and constant surface tension-is to be believed, then a contradiction ensues since at the triple point or contact line the velocity is both zero and non zero [1]. Attempts to solve this paradox and make progress on the issue abound [2,3,4]. One of the most popular is the assumption of Navier slip [5], but in general all solutions to the paradox amount to the introduction of a small length scale l µ below which the continuum model ceases to be valid a e-mail: ugis@mech.kth.se arXiv:2003.12246v1 [physics.flu-dyn] 27 Mar 2020 as discussed by Voinov [6]. Cox [7] extended Voinov's theory to arbitrary viscosity ratios and contact angles. Since then, many theoretical endeavours has been directed towards solving this problem [8,9,10,11,12,13,14,15,16,17] and the effort continues. Nevertheless, the microscopic scale physics remain somewhat mysterious as described for example in the review of [18]. Indeed, and beyond the paradox described above, there are many uncertain features of the nanoscopic flow and interface shape: there is uncertainty about the value of the contact angle at the smallest scale, and about the nature of the deviation from equilibrium, the effect of molecular forces, the presence of evaporation, etc. Despite recent advancements [19], experiments have difficulty providing interface shape and velocity field data for moving contact lines below the micron scale. From the applied mathematical point of view, Navier slip regularisation leads to an approximation in which the curvature still diverges logarithmically at the contact line and to a contradiction if the velocity field is continuous [20,21]. Moreover certain fluid and surface combinations have very small slip lengths, below the nanometer scale [22,23]. In such systems, if one considers the problem at smaller and smaller scales, other molecular effects will become relevant before the slippage effects.
Thus a full experimental characterisation at the nanometer scale is difficult and it very well may be that, at least for some time, only molecular dynamics (MD) "numerical experiments" will provide the insight necessary to understand which regularisation is adequate, and particularly so when the slip length is small and other effects than slip are required for a well-posed problem. However the largest possible MD simulations in three dimensions are limited in physical size to only tens of nanometers in each direction. This would make it impossible to perform asymptotic matching with numerical solutions of the Navier-Stokes equations. Indeed, in order to perform the matching over all scales necessary, the Navier-Stokes equations would have to be solved down to the nanometer scale. Moreover, for a 1 mm droplet, refining the calculation to the nanometer scale would require a range of scales of 10 6 , which implies the impossible-to-attain number of 10 18 grid points. Even restricting the computations to two dimensions of space, 10 12 grid points sit on the borderline of currently feasible computations (and not all problems warrant the use of a supercomputer). This creates the necessity either for investigations limited to much smaller scales, or for an intermediate-scale model, between the molecular scales and the scale of the sharp-interface model. Such an intermediate scale model may be provided by the diffuse interface approach, in which the hypothesis of a sharp transition between phases is replaced by that of a smooth transition over a small length scale . In the words of Len Pismen [24] "Multiscale methods employing different techniques at disparate length and time scales is the only feasible way to overcome the scale gap in practical computations". Simultaneous usage of MD, a diffuse interface model and a sharp interface model would be an implementation of this program, with the diffuse interface approach playing a key role at the intermediate scale. In such a framework, MD would help design and define parameters for the diffuse interface model, which in turn would perform a similar service for larger scale sharp interface models.
Beyond offering an intermediate scale, diffuse interface models open the possibility of regularising the contact line problem in an efficient and physically meaningful way. Indeed, together with MD, they have led to a generalisation of the Navier slip boundary condition (so called GNBC) that takes into account the "uncompensated Young's stress" 1 and makes the contact line motion problem well posed [26,27]. On the other hand, diffuse interface models [28,29] introduce at least two additional length Fig. 1. Illustration of the two-phase forced wetting configuration. Water parameters are determined from the molecular dynamics water model, at temperature T = 300 K. The vapour parameters for phase field and volume of fluid models are 1000 times smaller. The chosen parameters lead to capillary number Ca = 2µwU b /σ = 0.05. scales, the interface thickness = β/α and the diffusion length l d = √ µ l M (the reader is referred to section 3 and appendix B for a presentation of the diffusion-based phase field (PF) model parameters). The number of length scales is different, for example, for van der Waals model [30] or Cahn-Allen model [31]. Since the newly introduced length scales come in addition of the slip length l s there are potentially three length scales that can play an important role in the physics of the contact line. In addition a diffuse interface model has several other parameters specifically related to the contact line: the so called contact line friction [32,28], the surface energies (that lead to the equilibrium contact angle through Young's relation) and their spatial distribution. This makes the selection of the parameters rather difficult: there are three length scales and two or more other parameters to adjust or select. In a way, this is the price to pay to have a more realistic continuum description of the nanoscale. This motivates our choice of a physical situation where two simplifications occur: small slip and negligible evaporation. Indeed, for water on silica far from the critical point, slip is relatively small and evaporation is moderate. This situation may be reproduced using MD based on the SPC/E water model, which allows us to capture the strong hydrogen bonds between the water molecules and silica molecules on the surface [22]. Many possible application targeted configurations can be investigated, such as capillary driven flows [33] or forced wetting flows [18]. For the purpose of this work, we use a two-dimensional water drop enclosed by two moving walls, as shown in Fig. 1 at sufficiently small capillary number Ca = 2µ w U b /σ = 0.05 for existence of steady state configuration [34,35]. Here, U b is wall velocity, µ w is water viscosity and σ is surface tension. The choice of water in contact with a silica-line substrate makes the work reported here markedly different from previous MD investigations of the contact line dynamics in the same geometrical configuration [26,36,18]. The chose geometrical configuration is also known for its simplicity: it involves a symmetrical setting with identical solid walls providing us with identical two phase interfaces. Using this configuration we report on a first attempt to match MD, PF, and volume of fluid (VOF) simulations. By matching we mean finding the PF or VOF parameters that best reproduces the steady shape of the interface and the steady velocity field from the MD. This procedure is, as far as we know, original. Several other studies of the correspondence between the PF parameters and MD were performed [37] but not the direct attack on the contact line dynamics problem we suggest here. This paper is organised as follows. In section 2 we describe the MD simulations and the measurements, later used for the benchmarking of PF and VOF. Then, in section 3 we describe the PF model we use, the matching procedure to reproduce the MD results and show how does PF and MD compare. The VOF model and comparison between the VOF and MD is shown in section 4. Next, in section 5 we discuss the presented results and some open questions that remain. Finally, in section 6 we draw conclusions from this study.
2 Molecular dynamics simulations of water over silica-like substrate Molecular dynamics describes the system on a molecular level. The simulated system consists of water molecules, which interact with other water molecules as in the system as well as the substrate through a specified force field. Molecule positions and velocities are then integrated over time and the results sampled to obtain a representation of continuum variables.

Setup
To represent the chosen water-vapour geometry (Fig. 1), a two-dimensional shear system with water in vacuum is created by placing a water slab between two walls, as shown in Fig. 2(a). After the initial equilibration step, some of the water molecules move to the void and essentially form a very sparse water vapour. The walls consist of rigid SiO 2 molecules which are neutral electric quadrupoles. Water cannot easily slip over this substrate due to the electrostatic interactions between water molecules, which are dipoles, and the substrate quadrupoles. Water forms hydrogen bonds with this substrate, a form of bond that is transient but strong enough to prohibit slip [23]. The surface-water interaction is tuned to yield equilibrium contact angle θ 0 = 97.5 • .
Data is collected inside bins of size 0.25 × 0.25 nm 2 along x and y. These bins form a regular grid covering the entire system. The binning is visualised along the y axis in  . Average water molecule velocity, mass and temperature is sampled inside all bins over a period of 50 ps, after which the data is stored and the bins are reset for the next sampling period. We perform four simulations from different starting configurations. These configurations are created by generating the initial velocity field with different seeds for the random generator. Before the shear is applied each configuration is allowed to relax over a period of 100 ps. More details about the procedure are available in appendix A.

System parameter measurements
Parameters for the MD system are reported in Fig. 1. Bulk water density ρ w , viscosity µ w and surface tension σ are measured in separate simulations of pure water. Following Quan et al. [26] the water slip over the substrate is characterised through a friction parameter β f = µ w /l s , where l s is the corresponding Navier slip length. The friction parameter β f = 5.9 µ w nm −1 is measured in a Couette flow setup using the boundary layer of water molecules in direct contact with the wall molecules. From this measurement, the Navier slip length l s = 0.17 nm is extracted. Finally, the interface width is taken as the length over which the density goes from bulk to vapour density. This changes over a range of three bins, giving = 0.75 nm.

Results
Shear experiments consist of two stages. Once the shear is applied the water slab deforms until reaching its equilibrium (steady) state. For the transient state we measure the separation ∆x(t) between the upper and lower contact line positions, which begins at zero in the original slab geometry and increases to x 1 in the steady state. The contact line (and liquid interface) positions are identified as the bins with a density half that of the bulk density. The mean separation for the four simulations is shown in Fig. 3(a). The steady state is reached in 5.9 ns, after which the data has been used to generate the mean result. We obtain a mean displacement ∆x = 5.87 nm.
In the steady state the slabs are characterised using the full interface shape x(y) and the flow field. We mirror the right interface along both axes and overlay it on top of the left interface. Thus the lower end represents the receding contact line and the upper end the advancing. The average of all interfaces from MD in the steady regime is shown in Fig. 3(b). To gain more detailed insight into the agreement between MD and continuum models, we introduce parametric coordinate s as illustrated in Fig. 3(b). The interface angle θ (s) is shown in the inset of Fig. 3(a). Finally, we extract an averaged flow field from the MD runs, which we will use to discuss the agreement between MD, PF and VOF (Figs. 4, 5, 6 and 9).

Phase field model
We consider a 2D phase field (PF) model of a binary mixture to model two regions of different densities and viscosities. The water and vapour phases in Fig. 1 are assumed to be incompressible and the interface between the regions to be diffuse, i.e. that quantities vary smoothly over the interface.

Governing equations for the binary mixture
The phase field model introduces a phase variable C(x, y) ranging from 1 in the water phase to −1 in the vapour. The derivation of the governing equations for the phase variable can be found in [38,29]. The phase field variable is governed by a convectiondiffusion equation in a yet undetermined flow field u as The diffusive flux J d = −M ∇φ is proportional to the gradient of the chemical potential, defined as φ = βΨ (C) − α∇ 2 C, with proportionality coefficient M , which is called the phase field mobility. Due to the assumption of an incompressible flow, the convective flux takes the simple form J c = u · ∇C. In the chemical potential, we have two parameters α and β, which are related to the surface tension as σ = 2 √ 2αβ/3 and the characteristic thickness of the diffuse interface as = α/β. In addition, it contains the derivative of the standard double-well potential Ψ (C) = (C + 1) 2 (C − 1) 2 /4. The motion of the fluids is described by the incompressible Navier-Stokes equations with variable density and viscosity. We define the density and viscosity as respectively. Here ρ w and µ w are the density and the viscosity of the water, and ρ g and µ g are the density and the viscosity of the vapour component. The Navier-Stokes equations then become where the volume force f σ = −C ∇φ corresponds to the surface tension force and acts over the diffuse interface region. This form of the surface tension forcing is the so called potential form [39], which uses a reduced pressure.

Boundary conditions for the phase field model
The convection-diffusion equation (1) is a fourth-order partial differential equation and requires two boundary conditions. First, we impose a non-equilibrium wetting condition [29,26,28], where µ f is a contact line friction parameter, having the same units as bulk dynamic viscosity, and Γ is the triple contact line between the water, vapour and solid wall.
Here, θ 0 is the equilibrium contact angle and g (C) = 0.5−0.75C +0.25C 3 a switching function describing a smooth transition from water to vapour side. The unit normal vectorn is defined for the fluid domain and is directed into the surrounding surface.
If one sets µ f = 0, the dynamic contact angle is always enforced to the equilibrium angle θ 0 . Non-zero contact line friction allows the dynamic contact angle to evolve naturally as a function of contact line speed. The second boundary condition for the phase function is zero diffusive flux of chemical potential through the boundaries, i.e. ∇φ ·n = 0. The fluid momentum equations are subject to zero wall-normal velocity, u y = 0. For the tangential velocity component, we first consider the classical no-slip condition, which is equivalent to setting the fluid velocity near the moving wall to the velocity of the wall, u x = U w . In this situation, the only mechanism through which the contact line can move is the diffusion of the phase field variable. We also impose the Navier slip condition, u x = U w − l s ∂ y u xny , where l s is the Navier slip length. Interestingly, it has been recently shown using basic kinematic considerations that the Navier slip condition is not the appropriate way to regularise the contact line region [21]. The final boundary condition which we investigate for the PF model is the generalised Navier boundary condition (GNBC). This boundary condition (using the friction factor β f for the slip velocity [26]) takes the form where the second term is the uncompensated Young's stress. We assume a constant slip length l s over the whole solid surface. This leads to a friction coefficient β f , which varies in the same way as the liquid viscosity. The equations are implemented in a finite-element solve. For more details see appendix B.

Comparison with molecular dynamics
The PF parameters contain quantities uniquely determined from MD (ρ w , µ w , σ, θ 0 , ) and quantities that are less known or fully unknown (µ f , M , ρ g and µ g ). The vapour properties are not known due to the small system size, which renders any measurements from the MD outside of the liquid phase impractical. Therefore we have fixed both ρ g and µ g in the PF (and also later in the VOF) simulations as small as numerically feasible, ρ g = 10 −3 ρ w and µ g = 10 −3 µ w , see also Owing to the fact that we have two independent parameters to fit, we set out to find a set of parameters for PF simulations, which would agree with MD results up to the accuracy we see in the MD. The fitting procedure, later referred to fit I, is as follows: 1. Adjust the contact line friction (µ f ) individually for advancing and receding contact lines to match the MD dynamic contact angles. 2. Adjust the phase field mobility (M ) to match the MD drop displacement.
Results are shown in Fig. 4(a,b). There is a local error in the interface angle (Fig. 4b) near the advancing contact line. The drop displacement error in comparison to MD does not exceed 2%, which we consider a very good match. The parameters needed to arrive with the corresponding phase field results are listed in Tab. 1.  With red line, we illustrate the isoline of C = 0, which we use as the representation of the interface, and the blue background field show the variation of phase field variable C.
0.070 0.070 0.075 Table 1. Summary of the unknown phase field parameters obtained using fitting procedure I.
10 nm × 10 nm close-ups of MD streamlines near the advancing and receding contact lines are shown in Fig. 4(c,e), respectively. For comparison, in Fig. 4(d,f) we show 20 nm × 10 nm close-ups of streamlines from PF with GNBC. We observe that PF streamlines cross the interface and extend several nanometers into the vapour. This is in contrast to the MD results and therefore the PF predictions obtained through fitting procedure I provide un-physical flow fields. The PF flow with no-slip and Navier-slip conditions show slightly worse agreement and are not reported. The extent of which PF streamlines cross from water to vapour in the PF results can be characterised with the Pe number (appendix B), quantifying the relative importance between convection and diffusion. From the last row in Tab. 1 we observe that all fits have resulted in a Pe which is much smaller than unity: a high diffusion regime, in contrast to the convection dominated MD data. The high diffusion leads to the many streamlines crossing the PF interface. If one would use a much larger Pe number (a smaller mobility) but keep other parameters the same, the obtained steady state drop displacement ∆x would be largely overestimated.
To obtain a more physical flow field, we devise another fitting procedure, in which we relax the requirement on the dynamic contact angle. Fitting procedure II is defined as: 1. For no-slip simulation, we select µ f = µ w and then vary M to match the drop displacement observed in MD. 2. For Navier-slip and GNBC boundary conditions, we adjust µ f to match the contact angle in the no-slip simulation and then fit M to match the drop displacement in the MD.  Table 2. Summary of the unknown phase field parameters obtained using fitting procedure II.
We carry out these fits for all considered boundary conditions. By investigating the obtained interface angle distribution (Fig. 5,a) we see that overall agreement is good, while the local error near contact lines is increased. This difference is, however, not visible in the interface shape: practically the same agreement as presented in Fig. 4(a) is obtained. However, investigating the flow field near the advancing contact line reveals much better agreement with the MD results, although there is a small overshoot through the interface of streamlines originating from close to the wall (for fair measurement of the overshoot, we have identified a single streamline in all the simulations, which originates 0.25 nm away from the wall within the drop). The improved agreement with the MD (and smaller overshoot of the streamlines) is due to a smaller phase field mobility (or larger Pe number, see Tab. 2). The reduced contact line friction is the reason why it is possible to use smaller phase field mobility. Reducing contact line friction leads to smaller displacement of the drop, which consequently allows to reduce the phase field mobility to increase the drop displacement back to the amount measured in the MD. These results suggest that as the phase field mobility M is reduced (or Pe number is increased), the flow field near contact line approaches the one observed in MD. Therefore we devise the final fitting procedure III, which is defined as: 1. Set the contact line friction µ f = 0 (fixing the dynamic contact angle to the equilibrium one). This gives minimum friction at the contact line. 2. Fit the phase field mobility M to obtain the drop displacement observed in MD results.
The obtained Pe numbers using fitting procedure III are given in Tab. 3. The obtained results for PF are shown in Fig. 6 along with Navier-slip results from fitting procedures I and II. For the fitting procedure III it is meaningless to use GNBC, since the GNBC is equivalent to Navier-slip boundary condition for µ f = 0. In Fig. 6(a) we see the increasing local error of angle near the contact lines, which remains localised in a thin region near both walls. By investigating the streamlines (Fig. 6,b-e) we conclude that indeed as the PF mobility is reduced, the overshoot of the streamlines in PF is reduced and consequently the agreement with MD improves. 2.9 6.0 Table 3. Summary of the unknown phase field parameters obtained using fitting procedure III.

Volume-of-fluid model
The VOF model is known to be well-suited for solving interfacial flows. The interface between water and vapour in VOF-in contrast to the PF-is reconstructed in a sharp manner. Therefore the VOF model a good candidate and is typically used to solve two-phase flows in macroscopic systems. In this section we investigate how accurately the model can capture the behaviour of the nanoscopic droplet.

Governing equations for the two-phase flow with VOF model
In the VOF method, the fluid momentum is governed by exactly the same equations as in the PF method, namely, incompressible Navier-Stokes equations with variable density and viscosity (3)(4). The difference in methods lies in three aspects. First, the governing equation for the concentration function C, which within VOF method is typically called volume fraction, is a convection equation instead of the convection-diffusion equation (1) for the PF concentration. Second, the volume fraction C varies from 0 in the vapour phase to 1 in the fluid phase and consequently fluid density and viscosity becomes ρ (C) = C ρ w + (C − 1) ρ g and µ (C) = C µ w + (C − 1) µ g , respectively. The third difference lies in the fact that the surface tension force within the momentum equation (3) takes form where κ is the local curvature of the interface, δ s is a discrete Dirac distribution function used to spread the surface tension force to the fluid mesh andn is the normal of the two-phase interface. The curvature and interface normal are approximated using the Continuum-Surface-Force approach, which states that The accuracy of the surface tension term is directly dependent on the accuracy of the curvature calculation. The height-function methodology is a VOF-based technique for calculating interface normals and curvatures. The interested reader can find more details on the VOF method in appendix C.

Boundary conditions for the VOF model
For the VOF model used in this work, we use a similar set of boundary conditions as for the PF model. As we have two sets of unknowns (flow field and volume fraction), we need boundary conditions to determine both. Due to the fact that we are interested in steady regime only, we use a constant angle wetting condition. The angle is imposed through modification of height functions near the boundary and essentially sets the orientation of interface normal as illustrated in Fig. 7. The imposed contact angle affects the overall flow calculation in two ways. First, it defines the orientation of the VOF interface reconstruction in cells that contain the contact line and, second, it influences the calculation of the surface tension term by affecting the curvature computed in cells at and near the contact line. The second boundary condition is the velocity condition for the fluid flow. Although the VOF method does not allow for any diffusive transport of the contact line, there is always some mesh dependent numerical slip of the contact line [41]. This, however, does not provide accurate control over the results, therefore the explicit specification of a boundary condition compatible with a moving contact line is desired. As we have observed from PF result comparison with MD, the Navier-slip condition alone is sufficient to match the MD results, therefore we restrict the study on VOF only on the Navier-slip condition. Furthermore, since the observed slip length from the MD is very small -recall l s = 0.17 nm -we made choice to localise the slip condition in the near vicinity of the moving contact line and away from the contact line impose the wall velocity u x = U w . The localised Navier-slip condition takes form Here locality of the slip condition is enforced using the bell function The bell function takes as an argument the distance from the contact line d = x−x CL and the width of the bell function b . For coordinates further than b away from the contact line the function is set to zero to recover the no-slip condition. The wall normal velocity component is set to zero u y = 0, same as in PF. For more details on implementation of the velocity condition in the VOF model, see appendix C.

Comparison with molecular dynamics
To obtain results from the VOF method, we carry out a fitting procedure, described as follows.
1. Fix the width of the bell function b = 3.91 nm or five grid sizes in order to capture the variation of the slip length with sufficiently many points. 2. Fix the advancing and receding contact angle to θ a = 101 • and θ r = 93 • to represent the dynamic angle in the steady regime. 3. Adjust the magnitude of the local slip length l s to match the drop displacement from the MD.
Note that for the wetting condition, ideally we would like impose some relationship or governing law relating contact angle and contact line velocity. However, because we are looking only at steady regime currently, we bypass the implementation of a proper dynamic contact angle model for simplicity. Interf. shape  Using this procedure, we have found that the local slip length, which provides the best match between the VOF and the MD is l s = 8 nm. In Fig. 8 we show the interface shape and interface angle comparison between the VOF, PF and MD results. From Fig. 8(a) we observe that interface shapes are practically indistinguishable. The interface angles from VOF, however, show similar local errors as the PF angles. Interestingly, for the first and last interface angle we observe a rather large jump in the VOF result, which might be an artefact of the constant contact angle imposition on the first mesh cell next to the wall.

The role of the velocity boundary condition
As we have observed in the results of the PF simulations, the velocity boundary condition applied to the fluid momentum equation does not seem to be important for global measures, such as the drop displacement. Regardless of the imposed boundary condition, the drop displacement can be captured correctly. Therefore, for the VOF investigations, we have focused only on local Navier-slip boundary condition.
As for the interface shape, which is a more local measure, we have concluded that both benchmarked continuum models (PF and VOF) had their shortcomings and the obtained interface angle always has some local errors near contact line. Notably, the PF has been successful in describing the interface angle at the receding contact line (Fig. 4,b), but it always produced local error near the advancing contact line. However, when trying to match the interface angle as close to the MD as possible, we observe that the PF model provides non-physical flow field near the interface (Fig. 4,c-f), with streamlines crossing the interface and continuing in the vapour phase. This inaccuracy can be averted by allowing larger local errors for the interface angle (Figs. 5,a and 6,a), which then reduces the streamline crossing into the vapour phase significantly (Figs. 5,c-e and 6,c-e). Note that the local error can be only observed when plotting the interface angle (Figs. 5,a and 6,a), while the interface shape is globally indistinguishable (Fig. 4,a).
The improvement of the flow field was obtained irrespective of chosen boundary condition (Fig. 5,c-e), while the choice of the boundary condition introduce minor local improvements in either the flow field (if we demand the same accuracy of local interface angle near contact line) or in the local interface angle near contact line (if we demand the same accuracy of local flow field, set by the chosen PF mobility). The reason for improved flow field is that adding more slippage (first with Navier-slip, and then with GNBC) allows us to use smaller PF mobility M and consequently obtain streamlines that follow the two phase interface more accurately.
For the best possible representation of the flow field, we have set the dynamic contact angle to the equilibrium values, which allowed us to use the smallest PF mobility parameter. This choice renders the GNBC condition obsolete. By comparing no-slip and Navier-slip conditions, Navier-slip has produced the most accurate streamlines with still acceptable representation of the interface shape. However, it is also possible to use the no slip condition with minor decrease of the flow field accuracy. Therefore our results indicate that the GNBC condition (or even Navier-slip condition) is not necessary for physically acceptable predictions for the flow system considered in the present work using the PF model. In addition, increase of the diffusive transport near the contact line seem to always worsen the flow field prediction, which suggests that sharp interface method could bet the best possible choice for modelling the present system.
For the accurate representation of the drop displacement in the VOF, we have fixed the width of the local Navier-slip condition and adjusted the amplitude. We have observed that using one particular slip length locally near the contact line allowed us to reach exactly the same drop displacement as MD. However, similar as for PF, the interface shape has small local errors in angle (Fig. 8), while interface shapes are globally indistinguishable.

Wall location and slip condition
To gain further insight into which boundary condition would be the most appropriate to model the system consisting of a water drop between two no-slip plates, we look at velocity distribution as a function of a distance from the contact line. The MD results are obtained as an average over the 4 runs centred to the moving contact line. The MD results are shown in Fig. 9 with black crosses.  For comparison we show PF with Navier-slip condition, zero contact line friction µ f = 0 and Pe = 6, which yielded the best fit of streamlines. Sampling the velocity field along the wall (Fig. 9, red dashed) provides unsatisfactory agreement with the MD both in the near vicinity of the contact line and also deeper in the drop. However, if the sampling location is moved 0.75 nm away from the wall (Fig. 9, green dashdotted), very good agreement is observed between MD and PF. Similar observations we have also made with the VOF method, see Fig. 9, blue dotted line, in which we present VOF results sampled 0.39 nm away from the wall. Moving the sampling location up or down would provide similar influence as in the PF. Here, an interesting observation is that the VOF flow field seems to be less sharp compared to the PF, which probably stems from the fact that mesh resolution near the interface in VOF (∆s V OF = 0.782 nm) is much coarser than the resolution of the PF model (∆s P F = 0.195 nm, appendix B.2).
These observations raise two important questions about the presented benchmark. First, the Navier-slip condition-based on the slip length measured from the MD (l s = 0.17 nm)-does not seem to produce accurate slip velocity even more than 10 nm away from the contact line. This naturally leads to a question about where the solid wall in the continuum modelling viewpoint should be located, compared to the MD molecular picture (Fig. 2). The wall location is important for application of the Navier-slip condition, as displacement of the wall would cause a direct influence on the needed slip length. The velocity profile agreement at 0.75 nm distance from the wall in the continuum simulation seem to suggest that for the best agreement between MD and PF, the wall in the PF simulations should be shifted downwards by 0.75 nm corresponding to the center of the first bin in Fig. 2(a). The other option would be to use much larger slip length at the wall (larger by a factor of 4 compared to what is currently measured in the MD). The same discussion applies also to the VOF results. The slip and wall location therefore remain open questions.

Diffusion of phase field model
The Péclet number in theory exists in the MD system, because one can define a self-diffusion or a heat diffusion coefficient. For the currently used MD water the self diffusivity is [22] D w = 2.3 · 10 −9 m 2 /s.
However, the MD system has a single species and thus has no mass diffusion. By "mass diffusion" we mean that one species diffuses relative to the other, that is if you consider the mass m 1 of one of the species, or its density ρ 1 , it obeys a diffusion equation.
The PF model in the current work, on the other hand, is incompressible and consequently does not contain any state, energy or heat equation and does not model any self diffusion. The PF model is, however, based on the diffusion of the phase field variable C, see equation (1), and consequently models mass diffusion of liquid. This consequently leads to mismatch between the parameters used in the PF model and diffusion properties determined from the MD system. The Péclet number determined using the self diffusion of the water, relative wall velocity and height of the water drop is P e = 2U w L D w = 20 · 50 · 10 −9 6 · 2.3 · 10 −9 = 72.5, while maximum value we have been able to use in PF simulations and represent the MD results accurately is an order of magnitude smaller, i.e. P e = 6 (see Tab. 3). In our opinion, the reason we need to have a moderate P e number in the PF is twofold; (i) the P e number can not be tow small because we need to capture the phenomena of very small cross-interface mass transport (equivalent to evaporation), but on the other hand (ii) P e number can not be too large because it also controls the amount of diffusive transport very close to the contact line and controls the finally obtained drop displacement. Nevertheless, due to differences in physical effects MD and PF capture, a perfect agreement in results should not be expected.

Conclusions
We have carried out simulations of steady water drop sheared between two moving plates using of MD, PF and VOF methods. The MD simulations allow us to probe the detailed physical picture of the moving contact line and also global measures such as drop displacement amplitude, interface shape and flow field.
By comparing the results from PF simulations to MD we have observed that a perfect match is not possible. There is however a choice of parameters and boundary conditions which provide a reasonable approximation of MD. We observed a trade-off between accuracy of interface shape and accuracy of the flow field. The best match with respect to the flow field provided the largest error in the interface shape and vice versa. Particularly interesting conclusion from this study is that the exact boundary condition for the flow in the PF simulation does not seem to play a major role. On the other hand, we have observed that for accurate predictions there is an upper limit of the Pclet number one can employ (Pe = 6.0). This is an important insight if one considers using PF model as a substitute for nanoscale simulations, as proposed by Kronbichler and Kreiss [42]. This is also in contrast to previous drop spreading MD and PF comparisons [22], where Pe = 1400 was used. This suggests that the PF parameters analysed here are not universal and depend on whether the contact line is advanced through forced wetting or capillary spreading. We have also concluded, that the GNBC, which in the literature is proposed as a solution for the moving contact line problem, is not necessary to model the selected configuration.
Based on the results from PF and MD comparison, we have decided to only investigate a localised Navier-slip boundary condition for the VOF method. With this approach, we managed to capture the drop displacement accurately, while the interface shape had similar local errors as observed in the PF model.
By comparing the bottom water layer velocity from MD with the results of PF and VOF, we have identified an interesting open question about the accuracy of the slip condition and the wall location. It is possible that from the perspective of the continuum modelling of such small systems, the wall location might play even more important role than the chosen boundary condition. We have discussed the role of diffusion in PF model and the physical phenomena it models, which is different to what happens in MD simulations, therefore a perfect agreement should not be expected.
There are also aspects, which has not been investigated in the present work. Although the interface width has been set based on insight from MD results, in general more thorough look at the interface thickness influence on nanoscopic PF results could be interesting. Also, the width of the localised Navier-slip condition in the VOF model and its influence on the agreement between VOF and MD has not been investigated. Additional investigations would contain a comparison of the transient behaviour, not only the steady regime. This could be a very challenging problem, because if the steady regime is not perfectly described, then capturing the transient correctly would pose even harder restrictions. A Details of molecular dynamics simulations The system described in section 2 and figure 2 consists of two parts: the water slab and the moving walls. The water molecules forming the slab are modeled as SPC/E water [43]. This is a relatively cheap water model which retains important qualities of real water: the three-point structure and a dipole moment. These features gives the ability for the molecule to form hydrogen bonds with other molecules, which leads to a very low slip length (l s = 0.17 nm, which is smaller than a molecular diameter) and high surface tension.
The walls are built of rigid SiO 2 molecules. Partial charges q Si = −2q O are set to the atoms, making the molecules overall charge neutral but with a quadrupole moment. This electrostatic interaction works with the water molecules to create the aforementioned hydrogen bonding. The charge value q O = −0.40e is used to yield the static contact angle θ 0 = 97.5 • . The molecules are set in an fcc (111) structure with a spacing of 0.45 nm. They are kept in the desired formation by applying harmonic restraints with spring constant k = 10,000 kJ mol −1 nm −2 to the oxygen atoms of each molecule. The molecules are thus free to rotate in the surface normal plane. To move the walls at a desired speed we shift these restraining positions which pulls the molecules along with them while still allowing for natural thermal motion. Remaining details about the wall are described in [22].
Simulations are performed using Gromacs 2019 [44] in double precision. Atomic positions and velocities are updated using the leap-frog integrator with a time step of 2 fs. Non-bonded van der Waals interactions are treated fully up to a cutoff of 0.9 nm. Coulomb interactions are treated using PME electrostatics which interact over infinite range, including all periodic images of the system. Periodic boundary conditions are applied along x and z. Along the y axis reflecting walls are placed at each end to contain molecules to the system, although water molecules are kept from reaching these reflecting walls by the physical SiO 2 walls. A velocity rescaling thermostat is applied to the SiO 2 walls to dissipate excess energy with a time scale of 10 ps and keep the system at simulation temperature T = 300 K. Outside of the initial 100 ps equilibration period, the thermostat is not applied to the water molecules which can only dissipate energy by interacting with the walls through friction. The size of the full simulation domain along x is 150 nm, twice that of the water slab. Whereas a continuum model can be purely two-dimensional, a molecular system is naturally three-dimensional. A quasi-2D system is created with a system thickness of 4.667 nm in addition to the domain width and height. With this thickness, a width of 75 nm and height 50 nm, the water slab for the simulation is constructed with ∼ 580,000 water molecules.

B Details of phase field simulations
In this appendix, we describe in more details the non-dimensional governing equations as well as numerical implementation of these equations.

B.1 Dimensionless equations and dimensionless numbers
In order to render the dimensional governing equations (1,3,4) and the corresponding boundary conditions dimensionless, we choose the relative plate velocity U = 2 U w , the drop height L and the ration U/L as characteristic scales for velocity, length and time, respectively. Under this assumption, the convection-diffusion equation for phase field variable C (1) becomes where all variables are now dimensionless. This procedure gives rise to Péclet number and Cahn number, The Péclet number is quite important as it provides a measure of relative importance between convective and diffusive transport of the PF concentration function C. This number we use in the main text to report the needed PF mobility values for the agreement between PF and MD, see Tabs. 1, 2 and 3.
The incompressible Navier-Stokes equations (3)(4) in non-dimensional form becomes where again all variables now are non-dimensional. During this process, two new dimensionless numbers are introduced, namely Reynolds and capillary numbers, as The density and viscosity is normalised with respect to the parameters of water, and therefore in the dimensionless setting becomes where ρ g and µ g are the density and the viscosity of the vapour part, as shown in Fig. 1. Finally, the dimensionless wetting condition becomes Note that the numerical pre-factors 2 √ 2/3 appearing in certain places of the dimensional governing equation stems from relationship between the dimensional phase field constants α and β and the surface tension, i.e., σ = 2 √ 2αβ/3.

B.2 Numerical implementation
The dimensionless governing equations introduced in the appendix B.1 are linearised, cast into the weak form and solved using a symbolic finite-element toolbox femLego [45], which allows easy specification of finite-element weak form to solve. The solver is based on finite-element package deal.II and includes adaptive mesh refinement to resolve the sharp transition of phase field variables over a thin interface [46]. Linear elements were used for the phase field variables, while the fluid flow was resolved using Taylor-Hood elements (quadratic for velocity and linear for pressure). Mesh resolution was selected after a refinement study based on PF simulation with a no-slip condition. The final resolution selected and used for all other simulations is ∆s 1 = 3.125 nm far from the interface, and down to ∆s 2 = 0.195 nm within the interface region. Constant time step was used through the simulation as ∆t = 0.002 dimensionless time units.

C Volume of fluids
In this appendix, we describe in more details the numerical implementation of the solver, the height functions employed for the VOF model, as well as the implementation of velocity boundary condition.

C.1 Numerical implementation
In our study of the sheared droplet system, we used the free software Basilisk, successor of Gerris, developed at the Institut Jean le Rond d'Alembert (Sorbonne Universit) [40], [47], [41], [48]. The incompressible Navier-Stokes equations are discretised using the finite volume method and are solved using second order Bell-Colella-Glaz projection scheme [49]. The solver is coupled with VOF method for interface tracking. For obtaining the solution, we have used a uniform mesh spacing of ∆s = 0.782 nm. The height-function methodology is a VOF-based technique for calculating interface normals and curvatures. About each interface cell, fluid "heights" are calculated by summing fluid volume in the direction most normal to the interface. In 2D, a 7 × 3 stencil around an interface cell is constructed (Fig. 10) and the heights are evaluated by summing volume fractions horizontally, i.e.

C.2 Height functions
where ∆ is the grid spacing. The heights are then used to compute the the interface normaln and the curvature κ aŝ n = (∂ x h, −1) and κ = respectively. Here, ∂ x h and ∂ 2 xx h are discretised using second-order central differences.

C.3 Implementation of velocity boundary condition
In order to impose the velocity boundary condition at the wall in the selected numerical implementation, one has to work with discretisation stencils near the solid boundary. These stencils will extend beyond beyond the wall and make use of so called ghost points (Fig 11). The stencil values outside the domain (ghost values) need to be initialised. These values are set in order to provide the discrete equivalent of the NBC as where u x [ghost] is the tangential velocity at the ghost cell, u x [ ] is the tangential velocity of the cell inside the domain and ∆ is the grid spacing. This implementation is used to impose the locally specified Navier-slip condition, presented in equation (10).