Sub-Nanometer Interfacial Hydrodynamics: The Interplay of Interfacial Viscosity and Surface Friction

TL;DR

This paper develops a framework to accurately model interfacial flow at the nanoscale by incorporating position-dependent surface friction and viscosity profiles derived from molecular dynamics simulations, revealing exponential relationships with adhesion work.

Contribution

It introduces a novel approach to account for interfacial viscosity and surface friction in nanofluidic flow modeling, bridging molecular dynamics and continuum hydrodynamics.

Findings

Power-law relationships among friction, viscosity, and depletion length.

Exponential dependence of interfacial properties on work of adhesion.

Validated framework for sub-nanometer interfacial flow modeling.

Abstract

For an accurate description of nanofluidic systems, it is crucial to account for the transport properties of liquids at surfaces on sub-nanometer scales, where classical hydrodynamics fails due to the finite range of surface-liquid interactions and modifications of the local viscosity. We show how to account for both via generalized, position-dependent surface-friction and interfacial viscosity profiles, which enables the accurate description of interfacial flow on the nanoscale using the Stokes equation. Such profiles are extracted from non-equilibrium molecular dynamics simulations of water on polar, non-polar, fluorinated, and unfluorinated alkane and alcohol self-assembled monolayers of widely varying wetting characteristics. Power-law relationships among the Navier friction coefficient, interfacial viscosity excess, and depletion length are revealed, and these are each found to be exponential in the work of adhesion. Our framework forms the basis for describing sub-nanometer fluid flow at interfaces with implications for electrokinetics, biophysics, and nanofluidics.

Equations

$$F_f=-\lambda u_{\mathrm{slip}},$$

$$u_{\mathrm{slip}}=u(z_{\mathrm{int}})=b{\frac{du(z)}{dz}|}_{z=z_{\mathrm{int}}},$$

$$F_f=\int \mathrm{d}z{f}_f(z).$$

$$f_f(z)=-\ell (z)u(z).$$

$$\lambda u_{\mathrm{slip}}=\int \mathrm{d}z\ell (z)u(z).$$

$$\lambda =\int \mathrm{d}z\ell (z),$$

$$u_{\mathrm{slip}}=\frac{\int \mathrm{d}z\ell (z)u(z)}{\int \mathrm{d}z\ell (z)},$$

$$f_a(z)=\frac{F_a{\rho }_{\mathrm{liq}}(z)}{\int \mathrm{d}{z}^{\prime }{\rho }_{\mathrm{liq}}(z^{\prime })},$$

$$f(z)=f_a(z)+f_f(z).$$

$$f(z)=-{\partial }_z\eta (z){\partial }_{z}u(z),$$

$$\eta (z)=-\frac{1}{{\partial }_{z}u(z)}{\int }_{z_0}^z\mathrm{d}{z}^{\prime }f(z^{\prime }),$$

$$\ell (z)u(z)=f_a(z)+{\partial }_z\eta (z){\partial }_{z}u(z).$$

$$W={\gamma }_{lv}\left(1+k\right),$$

$$\lambda \propto e^{A(k+1)},$$

Full text

{tocentry}

The nanoscale flow of liquids at surfaces, especially water, is of crucial importance in electrokinetics 1, 2, 3, 4, as well as cell biology and biophysics 5, 6, such as for transport in transmembrane pores 7, bacterial motility 8, and hydration layers around biomolecules 9. It is the cornerstone of nanofluidics, which is of great interest recently 10, with direct applications including carbon nanotube technologies 11, 12, 13, 14, fabricated nanopores and nanoslits 15, 16, clean energy 17, 18, nanofiltration 19, 20, and biological/medical applications 21, such as the manipulation of DNA 22, intracellular delivery of biomolecules 23, and the analysis of cells, vesicles, and viruses 24. We seek to understand, down to the sub-nanometer scale, how liquids flow adjacent to surfaces, for which the relevant transport properties are the viscosity and surface–liquid friction. As one moves to the sub-nanometer scale, surface-to-volume ratios grow very large, which means that interfacial properties (e.g., surface–liquid friction), become as important as bulk properties (e.g., viscosity) 25, 26. The surface–liquid friction is characterized for steady-state flows by the Navier friction coefficient, $\lambda$, defined via

[TABLE]

where $F_{f}$ is the surface–liquid friction stress, and $u_{\rm slip}$ the slip velocity of the liquid, i.e., the relative velocity of the liquid directly adjacent to the surface 27. Alternatively, the slip may be cast in terms of the slip length $b$, which satisfies $\lambda=\eta/b$, where $\eta$ is the shear viscosity of the liquid 28. In the macroscopic context, the slip length is useful for describing low-friction interfaces, such as gases 29, or very-high-viscosity liquids, such as polymer melts30, at solid surfaces. For liquids like water, slip lengths are typically on the order of nanometers, and the effect of slippage on flow increases with $b/h$, the ratio of the slip length to the channel size 25, so slip lengths for macroscopic flows are vanishingly insignificant, and a no-slip condition is sufficient 31. Despite this, small slippage at solid–liquid interfaces was measured as early as the 1950s in glass capillary channels 32, 33. More recently, liquids have been shown to flow through carbon nanotubes at rates multiple orders of magnitude larger than would a no-slip flow 11, 12, 13, 14. Note that eq 1 describes the linear-response regime, i.e. where the slip velocity is low; when it becomes large enough, the relationship is nonlinear 34, 35.

The Navier friction coefficient $\lambda$ and/or slip length $b$ can be extracted from non-equilibrium molecular dynamics (NEMD) simulations where shear is induced by applying external driving forces 34, 36, 37, 38, 39. A typical approach is to induce a steady-state Couette or Poiseuille flow and linearly extrapolate the tangential liquid velocity profile $u(z)$ past the interfacial position $z_{\rm int}$ to find the slip length $b$ via

[TABLE]

with $z_{\rm int}$ defined for example, as the (fixed) position of the solid surface 37 or as the Gibbs dividing surface of the liquid 40, 39. While this method works reasonably well for weak surface–liquid interactions and low-friction surfaces, problems arise for high-friction interfaces where $u_{\rm slip}$ and $b$ are small, difficult to measure, and sensitively dependent on the definition of the interfacial position, $z_{\rm int}$.

In another approach, equilibrium molecular dynamics (MD) simulations of a solid surface with an adsorbed liquid are carried out, the friction coefficient of the center-of-mass coordinate of the liquid or surface are extracted using Green–Kubo relations, and the Navier friction coefficient is found using Stokes’ equation, where a position- and time-independent viscosity $\eta$ is assumed 41, 42.

These models crucially ignore the non-locality of intermolecular interactions, which leads to two modifications of classical hydrodynamics: firstly, surface–liquid friction acts over a finite distance and is therefore delocalized along $z$, as opposed to acting just at a specific position $z_{\rm int}$, and secondly, the structure and behavior of liquid near an interface is modified, which means that shear viscosity near interfaces is in general position dependent. Indeed, viscosity near interfaces has been found to be enhanced when liquid interactions with the adjacent phase are strong, and reduced when they are weak 43, 40, 44, 45. Modified interfacial viscosity is of particular importance for nanoconfined liquids 46, 47. Thus, to fully describe interfacial hydrodynamics at the nanoscale, we take into account the position dependence of the surface–liquid friction coefficient $\ell(z)$ and interfacial viscosity $\eta(z)$, extracting both from simulations.

To this end, we carry out driven-flow NEMD simulations of liquid water slabs on self-assembled monolayer (SAM) surfaces with widely varying wetting and friction properties. Above each water slab is vacuum in which a water vapor phase may form. This open system geometry circumvents the strong sensitivity of viscosity to nanochannel width arising from the incommensurability of discrete liquid-molecule layers—an important phenomenon, but not the focus of this study 48, 49. The SAMs are comprised of decane (H-SAM), decane with the top eight carbons perfluorinated (F-SAM), or decanol ($\alpha$-SAM), where the partial charges on the OH groups are scaled by a factor $\alpha\in\{0,\,0.5,\,0.6,\,0.7,\,0.8,\,0.9,\,1\}$. The systems are illustrated in Figure 1(a), though the production simulation systems are larger than those shown.

The simulations are carried out in GPU-enabled single-precision GROMACS ${2023.3}$ 50, 51, using the leap-frog integrator 52 with a $2$-fs time step. The velocity-rescaling (CSVR) thermostat is applied to all atoms with a target temperature of $300$ K 53. Because longer force cutoffs are superior for modeling interfacial properties, the Lennard–Jones forces are modeled with force-switching between $1.9$ and $2$ nm 54, 55. Electrostatic forces are modeled using particle-mesh Ewald beyond a real-space cutoff of $2$ nm 56. SAM molecules are modeled using the OPLS All-Atom (OPLS-AA) force field 57, 58, 59 with selected dihedrals optimized 54. Water is modeled using the SPC/E water model 60. Further simulation details can be found in the SI Section S1.

The water is driven along the $x$-direction, tangential to the SAM, by a constant gravity-like force acting on all liquid atoms, i.e., the force on each atom is proportional to its mass. In the bulk, where the liquid density is constant, this results in a constant force density, and a quadratic flow profile with the vertex of the parabola at the liquid–vapor interface, i.e., a half-Poiseuille flow, as illustrated in Figure 1(b).

Figures 1(c–j) illustrate the workflow for extraction of surface-friction and interfacial-viscosity profiles from the driven-flow simulations, and resulting flows for the low-friction, fluorinated F-SAM (left column) and the high-friction, polar ($\alpha$=1)-SAM (right column). Figures 1(c, d) show snapshots of the simulated systems, zoomed in on the respective interfacial regions of interest. Figures 1(e, f) plot the extracted profiles $u(z)$, the liquid velocity, and $f_{f}(z)$, the friction force per unit volume on the centers of mass of liquid molecules due just to the surface (and not to adjacent liquid), as a function of the $z$-position of liquid molecule centers of mass. The slip at the low-friction, hydrophobic F-SAM is apparent in $u(z)$, i.e., $u(z)$ is clearly positive valued even approaching the surface. For the $\alpha$-SAM on the other hand, slip is not apparent, and $u(z)$ is convex near the interface, making the definition of a slip velocity or slip length difficult. The non-locality of $f_{f}(z)$ is apparent in Figures 1(e, f), with the surface–liquid friction force acting on the liquid over several Angstroms along $z$ for both surfaces. This illustrates why the surface–liquid friction force should not be ignored in viscosity profile calculations. The total surface–liquid friction stress $F_{f}$ is given by

[TABLE]

In the linear-friction regime, a local, position-dependent friction coefficient $\ell(z)$ is defined via

[TABLE]

Combining eqs 1, 3 and 4 yields

[TABLE]

Equation 5 holds for arbitrary profiles $u(z)$, including constant $u(z)$, with $u=u_{\rm slip}$, from which follows by comparison with eq 1 an expression for the Navier friction coefficient based on the microscopically defined friction profile,

[TABLE]

which, crucially, does not rely on the definition of the interfacial position $z_{\rm int}$. Equations 5 and 6 can be combined to give the corresponding slip velocity,

[TABLE]

which takes the form of a weighted mean over $u(z)$. Thus, from eq 4, it follows that the Navier friction law, eq 1, holds for $\lambda$ and $u_{\rm slip}$ as defined in eqs 6 and 7 (see the SI Section S2 for a more complete derivation). Surface-friction profiles $\ell(z)$, calculated via eq 4, are plotted for the F-SAM/water and ($\alpha$=1)-SAM/water systems in Figures 1(g, h). In order that the reader is able to discern the small differences in shape and position between $f_{f}(z)$ and $\ell(z)$, the forces $f_{f}(z)$ are also reproduced in Figures 1(g, h) as translucent orange curves with an arbitrary scaling factor. As one might expect, the friction profile for the very hydrophilic ($\alpha$=1)-SAM is (two) orders of magnitude higher than that of the hydrophobic F-SAM.

As mentioned above, viscosity is modified at interfaces. Let $f_{a}(z)$ denote the applied driving force density acting on the liquid in the $x$-direction,

[TABLE]

where $\rho_{\rm liq}(z)$ is the liquid density profile and $F_{a}$ is the total applied driving stress. Then, the total external force density on the liquid, $f(z)$, can be decomposed as

[TABLE]

Water on the nanoscopic scale has a low Reynolds number, so its motion is described well by the linear Stokes equation, 27, 61, 62

[TABLE]

where a local, position-dependent shear viscosity $\eta(z)$ is assumed (see the SI Section S3 for derivations of the Navier–Stokes and Stokes equations for position-dependent viscosity). Letting $z=z_{0}$ denote a position below the liquid phase, integration of eq 10 from $z_{0}$ to an arbitrary $z$ gives

[TABLE]

a more thorough derivation and discussion of which can be found in the SI Section S4. Viscosity profiles $\eta(z)$, calculated from eq 11, are plotted for the F-SAM/water and $\alpha$-SAM/water systems in Figures 1(g, h) (see the SI Sections S5 and S6 for details and plots on extraction of $\ell(z)$ and $\eta(z)$). Each viscosity profile converges to a bulk value $\eta_{b}$, which agree well between the two systems.

Several previous publications take a similar approach to eq 11, but take $f(z)\propto\rho_{\rm liq}(z)$, which, crucially, does not account for the surface–liquid friction stress 63, 64. Others take very different approaches, including measuring liquid shear between two oppositely moving surfaces 43, 65, 40 or in flow driven by electric fields 66, via the local pressure tensor 67, 68, 69, or via the shear and momentum flux 70, none of which, however, correctly separate surface–liquid stress and viscous effects. Ref. 71, on the other hand, calculates interfacial viscosity profiles for Lennard–Jones fluids via pairwise fluid-fluid interactions, which circumvents surface friction contributions.

Substituting eqs 4 and 9 into the Stokes equation, eq 10, gives

[TABLE]

Having extracted $\ell(z)$ and $\eta(z)$ for a specific surface–liquid pair at a given temperature and pressure, $u(z)$ may be found by solving eq 12 numerically for an arbitrarily chosen applied force profile $f_{a}(z)$ (see the SI Section S7). Figures 1(i, j) compare velocity profiles extracted from driven-flow simulations with those calculated by solving eq 12. For the liquid–vapor boundary, it is assumed that $\partial_{z}u(z)=0$, while at the surface, it is assumed that $u(z)=0$. The overall modeling of the flow is very accurate, which validates our model, including the locality of $\ell(z)$ and $\eta(z)$, and indicates that the system is in the linear-response regime for both the viscosity and friction (linear-response is also verified in the SI Section S8). The same comparison between the $\ell$–$\eta$ model and driven-flow simulation data is carried out for all systems studied in this work in the SI Section S9. In the SI Section S10, we show that failing to account for the position dependence of either $\ell(z)$ or $\eta(z)$ leads to models that fail to accurately capture the behavior near the interface.

The surface friction and viscosity profiles for all systems are collected in Figure 2. Wetting is a natural way to classify surface–liquid interactions, so we extract water contact angles of all systems from density profiles of equilibrium MD simulations of cylindrical droplets55, and extrapolate the droplet size to the macroscopic limit72, 73, 74, 54. The contact angle extraction is detailed in the SI Section S11. Figure 2(a) shows a simulation snapshot of such a cylindrical water droplet on an F-SAM, with the contact angle indicated schematically, and the resulting contact angles for all systems, which range from the very hydrophobic F-SAM ($\theta=125^{\circ}$) to the fully wetted ($\alpha$=1)-SAM.

Figure 2(b) shows the surface-friction profiles $\ell(z)$, with the corresponding mass density profiles shown above as a positional reference. As the friction varies by more than an order of magnitude, the curves for the low-friction surfaces are difficult to discern in the main plot, and the same data is plotted over a smaller range of $y$-values in the inset. The surface-friction profiles are delocalized over several Angstroms, with FWHM ranging from 1.1 to 1.8 Å across all systems. This reveals that, on sub-nanometer scales, friction should not be treated as acting only at a single position $z_{\rm int}$.

The viscosity profiles $\eta(z)$ are shown in Figure 2(c). Note that these all decay to zero in the region of the SAM, which has been observed before for Lennard–Jones fluids 71. If surface friction is ignored in eq 11, i.e., all flow behavior is attributed to the position-dependent viscosity alone, the resulting effective viscosity will typically be found to diverge near the interface (see the SI Section S12 for more on effective viscosities, including plots for all systems studied in this work). The viscosity profiles for all systems converge very near the same value in the bulk, as expected. For each system, the viscosity value in the bulk, $\eta_{b}$, is calculated as the mean over the profile in a window between 1.5 and 2.5 nm from the top carbon atoms. These are shown as dotted horizontal lines and they agree well, all falling within the range $\overline{\eta}_{b}\pm 0.007$ mPa s, where $\overline{\eta}_{b}=0.706$ mPa s is their mean value, which agrees well with the bulk value of $\eta_{\rm eq}=0.698$ mPa s, calculated via the Green–Kubo relation from an equilibrium MD simulation where the same water force field is used (see the SI Section S13 for details) 75, 76, 77, 78.

Alongside the viscosity profiles $\eta(z)$ in Figure 2(c), the viscosity dividing surfaces $z_{\eta}$ (where the viscosity excess vanishes) are shown as vertical dashed lines, and differ from the Gibbs dividing surface $z_{\rm GDS}$ by up to about an Angstrom in either direction. This reveals that $z_{\rm GDS}$ is not the relevant interface for liquid viscosity and highlights the importance of treating the position-dependence of the interfacial viscosity on sub-nanometer scales. We define the interfacial viscosity excess distance as $d_{\eta}=z_{\rm GDS}-z_{\eta}$. Thus, a positive $d_{\eta}$ indicates a viscosity excess at the interface, while a negative $d_{\eta}$ indicates a deficit. The viscosity excess grows significantly as the surfaces become more hydrophilic, with a shoulder forming in the profile at the first hydration layer nearest the surface, which likely results from the conformational rigidity of water hydrogen-bonded to, or otherwise strongly interacting with, the surface, preventing other water molecules from easily flowing past. One interesting aspect visible in all the systems is the apparent disagreement between the viscosity and density profiles: naively one might expect viscosity to increase with density, but we observe the main peak of each viscosity profile to be shifted bulkward relative to that of the density profile. This is likely related to water orientation and especially hydrogen bonding near the interface. Indeed, we show that there is a deficit of water-water hydrogen-bonded molecules near the interface for all surfaces in the SI Section S14.

In Figure 3, we explore the relationships among the Navier friction coefficient $\lambda$, the interfacial viscosity excess distance $d_{\eta}$, the wetting coefficient $k=\cos\theta$, and the depletion length $\delta$. The depletion length is the distance between the adjacent Gibbs dividing surfaces of the SAM and the water phases, as illustrated in Figure 3(a), where the Gibbs dividing surface of a phase is the thermodynamically relevant interfacial position (see the SI Section S15). For the SAMs, which are chemically inhomogeneous near the interface, the Gibbs dividing surfaces are calculated as the position where the excess of the packing density $\phi(z)$ vanishes (see the SI Section S16 for further discussion).

For partial wetting, the wetting coefficient $k$ is related to the areal work of adhesion of the surface $W$ via the Young-Dupré equation,

[TABLE]

where $\gamma_{lv}$ is the liquid–vapor interfacial tension of the liquid. Thus, we plot $k+1$ rather than $\theta$, and because eq 13 only holds for partial wetting, we do not include the data for the ($\alpha$=1)-SAM in plots involving $k+1$.

In Figure 3(b), $d_{\eta}$ is plotted over $k+1$ alongside a least-squares fit of an exponential function plus a constant, which we call $d_{\eta}^{0}$. From the fitted function we calculate $d_{\eta}(k=-1)=-1.25$ Å, i.e., there is an interfacial viscosity deficit of 1.25 Å in the dewetting limit, where $\theta=180^{\circ}$, which should be characteristic of the water model alone. Figures 3(c–e) are log-linear plots of $\lambda$, $d_{\eta}-d_{\eta}^{0}$, and $\delta$ over $k+1$, alongside exponential fits. The fits in Figures 3(c, e) are performed in logarithmic space.

From transition state theory 79, 80, it can be shown that $\lambda\propto e^{\beta U_{0}}$ for small driving forces, where $U_{0}$ is the barrier height of the corrugated energy landscape seen by a liquid molecule at a surface and $\beta=(k_{B}T)^{-1}$ is the inverse thermal energy. From there, it can be shown under certain approximations that

[TABLE]

where $A$ is a positive constant. Indeed, we find $\lambda$ to be exponential $k+1$ in Figure 3(c), where the fit gives $A\approx 2.3$. Previous works have derived other friction/barrier-height relationships, including $\lambda\propto e^{\beta U_{0}}/U_{0}$ and $\lambda\propto 1+\left(\beta U_{0}\right)^{2}/16$ 81, 43, 82. Fits of these functions to our $\lambda(k)$ data, a more detailed derivation of eq 14, and related discussion can be found in the SI Section S17.

Figure 3(d) shows the same fit as (b), but shifted by $-d_{\eta}^{0}$ and in log-linear, where the data indeed appear to fall on a straight line, indicating exponential behavior. Figure 3(e) plots the depletion length, and these data seem to follow an exponential decay.

Figures 3(f–h) examine the relationships among $\lambda$, $d_{\eta}-d_{\eta}^{0}$, and $\delta$. As these all appeared exponential in $k+1$, we expect them to relate to one another via power laws, and accordingly plot them in log-log. The data are fitted by taking the logarithm of both datasets and carrying out a linear fit. From Figure 3(f) we find $\lambda\propto\delta^{-2.9}$, from (g), $d_{\eta}-d_{\eta}^{0}\propto\delta^{-1.9}$, and from (h), $\lambda\propto(d_{\eta}-d_{\eta}^{0})^{1.6}$.

Thus, we have shown that both surface-friction and inhomogeneous interfacial viscosity profiles are important for understanding flow, and that properly decoupling the two allows for the calculation of useful quantities, such as the surface viscosity excess and Navier friction coefficient, without a priori recourse to arbitrarily defined interfacial positions. Going forward, these methods for accurately decoupling surface–liquid friction and interfacial viscosity will lend themselves well to the analysis of flow on sub-nanometer scales for a wide variety of systems, such as systems with liquid mixtures, nanochannels, and charged surfaces and/or electrolyte solutions, which are relevant for electrokinetics. The constitutive relations we find between friction coefficient $\lambda$, viscosity excess $d_{\eta}$, depletion length $\delta$, and wetting coefficient $k$ allow for the parameter-free modeling of interfacial flow at smooth surfaces of arbitrary polarity.

{acknowledgement}

The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for funding via Project-ID 387284271 (SFB 1349).

{suppinfo}

S1: Simulation Details. S2: Position-Dependent Surface–Liquid Friction. S3: Navier–Stokes and Stokes Equations for Position-Dependent Viscosity. S4: Calculation of Position-Dependent Viscosity. S5: Fits of Velocity Profile Tails. S6: Calculation of Friction and Viscosity Profiles for All Systems. S7: Numerically Solving the Stokes Equation. S8: Verifying Linear-Response Regime by Velocity Profile Fits. S9: Modeling Flow for All Systems. S10: Other Approaches to Modeling Flow. S11: Contact Angles. S12: Effective Viscosity Profiles. S13: Bulk Shear Viscosity from the Green–Kubo Relation. S14: Hydrogen Bonding Near the Interface. S15: Gibbs Dividing Surface. S16: Depletion Length. S17: The Friction–Wettability Relationship.