# A consistent treatment of dynamic contact angles in the sharp-interface framework with the generalized Navier boundary condition

Tomas Fullana, Yash Kulkarni, Mathis Fricke, Stéphane Popinet, Shahriar Afkhami, Dieter Bothe, Stéphane Zaleski

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Tomas Fullana, Yash Kulkarni, Mathis Fricke, Stéphane Popinet, Shahriar Afkhami, et al.. A consistent treatment of dynamic contact angles in the sharp-interface framework with the generalized Navier boundary condition. 2025. hal-05428626

HAL Id: hal-05428626

https://hal.science/hal-05428626v1

Preprint submitted on 22 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

# A consistent treatment of dynamic contact angles in the sharp-interface framework with the generalized Navier boundary condition

Tomas Fullana¹,², Yash Kulkarni¹, Mathis Fricke³, Stéphane Popinet¹, Shahriar Afkhami⁴, Dieter Bothe³, and Stéphane Zaleski¹,⁵

¹Sorbonne Université and CNRS, Institut Jean Le Rond d'Alembert UMR 7190, F-75005 Paris, France

²Laboratory of Fluid Mechanics and Instabilities, EPFL, Lausanne CH-1015, Switzerland

³Department of Mathematics, TU Darmstadt, Schlossgartenstraße 7, 64289 Darmstadt, Germany

⁴Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ, USA 07102

⁵Institut Universitaire de France, Paris, France

(Received xx; revised xx; accepted xx)

In this work, we revisit the Generalized Navier Boundary condition (GNBC) introduced by Qian et al. in the sharp interface Volume-of-Fluid context. We replace the singular uncompensated Young stress by a smooth function with a characteristic width $\epsilon > 0$ that is understood as a physical parameter of the model. Therefore, we call the model the “Contact Region GNBC” (CR-GNBC). We show that the model is consistent with the fundamental kinematics of the contact angle transport described by Fricke, Köhne and Bothe. We implement the model in the geometrical Volume-of-Fluid solver Basilisk using a “free angle” approach. This means that the dynamic contact angle is not prescribed but reconstructed from the interface geometry and subsequently applied as an input parameter to compute the uncompensated Young stress. We couple this approach to the two-phase Navier Stokes solver and study the withdrawing tape problem with a receding contact line. It is shown that the model allows for grid-independent solutions and leads to a full regularization of the singularity at the moving contact line, which is in accordance with the thin-film equation subject to this boundary condition. In particular, it is shown that the curvature at the moving contact line is finite and mesh converging. As predicted by the fundamental kinematics, the parallel shear stress component vanishes at the moving contact line for quasi-stationary states (i.e. for $\dot{\theta}_d = 0$), and the dynamic contact angle is determined by a balance between the uncompensated Young stress and an effective contact line friction. Furthermore, a non-linear generalization of the model is proposed, which aims at reproducing the Molecular Kinetic Theory of Blake and Haynes for quasi-stationary states.

**Key words:** dynamic contact line, Volume-of-Fluid method, Generalized Navier Boundary Condition, withdrawing plate, forced dewetting
# 1. Introduction

The phenomenon of dynamic wetting / dewetting requires a relative motion of the contact line, i.e. the triple line at which the liquid-fluid interface and the solid support's surface meet, against the solid wall. This fundamental process can be modeled in various ways. If the fluid interface and the contact line are modeled as a material surface and a material line, respectively, it is clear that the classical no-slip condition is incompatible with the dynamic wetting process. Mathematically, it has been shown that, for a material interface and contact line, the no-slip boundary condition leads to a discontinuity in the velocity as the contact line is approached. Because of that, a viscous fluid develops a non-integrable singularity at the moving contact line. This has been first shown in the seminal paper by Huh & Scriven (1971). Since then, various mathematical models have been developed to resolve the paradox in the continuum mechanical description. In the framework of diffuse interface models, where both the fluid interface and contact line have a finite width characterized by a smooth but rapidly varying order parameter, a motion of the contact line can be achieved by pure diffusion of the order parameter; see Jacqmin (2000). In this case, the motion is driven by gradients of the chemical potential and the contact line is not a material line with respect to the fluid particles. The Interface Formation Model due to Shikhmurzaev (1993, 2008) describes the dynamic wetting process using mass transfer between the bulk phases of the liquid and the interfaces between fluid and solid and fluid and gas. Hence, in this case, the contact line can move without hydrodynamic slip, displaying a rolling motion.

A commonly used approach in the sharp interface framework is to model the interface and the contact line as material objects and to allow for slip between the bulk fluid and the solid wall. The Navier slip condition states that the amount of tangential slip is determined by a balance between the tangential component of the viscous stress, described by the viscous stress tensor S, and a sliding friction force between fluid particles and the solid surface according to

$$ -\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) = (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel}, \quad (1.1) $$

where $\mathbf{n}_{\partial\Omega}$ denotes the unit outer normal to the solid boundary. This boundary condition introduces the slip length $\lambda := \eta/\beta$ as the key parameter, where $\eta$ denotes the viscosity of the liquid and $\beta > 0$ is a coefficient describing the (sliding) friction between the liquid molecules and the solid surface. Within the Navier slip model, the slip length can be interpreted geometrically as the distance below the solid surface where the linearly extrapolated tangential velocity vanishes. It is well-known that the singularity at the moving contact line is only partially relaxed by the Navier slip condition. A logarithmic divergence as a function of the distance to the contact line still exists for the curvature and the pressure, as pointed out by Huh & Mason (1977). However, the singularity is transformed into an integrable one and, hence, physically meaningful solutions are possible, at least for the macroscopic flow. The physical implications of the pressure singularity are debated in the literature. Shikhmurzaev (2006) argues that the pressure should remain finite because otherwise the model of an incompressible fluid would no longer be valid. On the other hand, it has been demonstrated that the Navier slip model is able to describe various wetting experiments in a satisfactory manner, see Fullana et al. (2020).

Besides the mobility of the contact line, the wettability of the solid surface is another key parameter for the physical system. It is usually characterized by the equilibrium contact angle $\theta_e$ that the fluid interface forms with the solid boundary in equilibrium. It can be computed from the surface tension of the liquid-gas, liquid-solid and solid-gas interfaces, using the
equation introduced by Young (1805), viz.

$$ \sigma \cos \theta_e + \sigma_{ls} - \sigma_{sg} = 0. \qquad (1.2) $$

While the latter equation can be easily deduced from variational principles, the dynamics of the contact angle is a much more complex problem and a large variety of empirical models exist. Notably, there is one fundamental relation for the dynamics of the microscopic contact angle $\theta_d$ in the limit of slow velocities of the contact line, which is shared by many of these models:

$$ -\zeta U_{cl} = \sigma(\cos \theta_d - \cos \theta_e). \qquad (1.3) $$

Here $U_{cl}$ denotes the normal speed of the contact line relative to the solid surface (positive for advancing and negative for a receding contact line) and $\zeta$ is the so-called “contact line friction” parameter. Equation (1.3) arises, for example, from the molecular kinetic theory of wetting in the limit of small capillary numbers, i.e. for a slow motion of the contact line (see Blake & Haynes (1969); Blake et al. (2015)).

Recently, Fricke et al. (2019, 2018) studied the fundamental kinematics of the contact angle transport and showed that the rate-of-change of the contact angle is fully determined by a directional derivative of the velocity field $v$ at the contact line, viz.

$$ \dot{\theta}_d = (\partial_{\tau} \mathbf{v}) \cdot \mathbf{n}_{\Sigma}. \qquad (1.4) $$

Here $\mathbf{n}_{\Sigma}$ denotes the interface normal vector and $\tau$ is the unit vector tangential to the interface and normal to the contact line (see Section 2 for more details). Notably, when applied to the full two-phase flow problem (assuming sufficient regularity of the solution) Equation (1.4) can be used to show that $\dot{\theta}_d$ is proportional to the tangential component of the shear stress that also appears in the Navier slip condition (see (2.16) in Section 2). In particular, this shear stress component vanishes if the contact angle does not change in time, i.e. in a quasi-steady state. In other words, (1.4) predicts an “apparent perfect slip” at the moving contact line if $\dot{\theta}_d = 0$. Indeed, indications of a vanishing shear stress in the vicinity of the contact line have been observed in molecular dynamics (MD) simulations by Thompson & Robbins (1989) and others. We will see that perfect slip in the sense of vanishing shear stress is possible within the GNBC model, which makes the model consistent with equation (1.4). On the other hand, Fricke et al. (2019) and Fricke & Bothe (2020) showed that the Navier slip model (1.1) with a contact angle boundary condition like (1.3) is not consistent with (1.4) and regular solutions, if they exist, show an unphysical behavior.

The “Generalized Navier Boundary Condition” (GNBC) was first described by Qian et al. (2003, 2006a,b) in the context of diffuse interface models and molecular dynamics. The key idea of the GNBC is to introduce the *uncompensated Young stress* as an additional force density into the constitutive relation (1.1). So, in this model, the slip velocity relative to the solid surface is a result of a balance between a sliding friction force, the viscous stress and the uncompensated Young stress. In a sharp interface and sharp contact line formulation, the GNBC can be written as (see Ren & E (2007), Gerbeau & Lelièvre (2009))

$$ -\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) = (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos \theta_d - \cos \theta_e) \mathbf{n}_{\Gamma}\delta_{\Gamma} \quad \text{on } \partial\Omega. \qquad (1.5) $$

Notably, the contact line delta distribution $\delta_{\Gamma}$ appears because, in the sharp contact line formulation, the Young stress is concentrated on a mathematical curve. Hence, Equation (1.5) should mathematically be understood as an equality of distributions. This delta function GNBC formulation is applicable in weak formulations of the two-phase flow problem where the contact line delta distribution will translate into an integral over the contact line in the
weak formulation (see, e.g, Gerbeau & Lelièvre (2009); Fumagalli et al. (2018)). On the other hand, there is no contact line delta distribution in the Phase Field formulation of the GNBC due to Qian et al., because the thickness of the interface and the contact line is a finite, physical model parameter in this case. Yamamoto et al. (2013, 2014) implemented the GNBC approach into a front-tracking-method and studied the dynamics of capillary rise in a tube. In this method, the contact line is transported by advecting the Lagrangian marker points without a prescribed contact angle. Then, the dynamic contact angle is evaluated and used to compute the uncompensated Young stress, which determines the slip velocity profile. Yamamoto et al. noticed that the viscous stress becomes negligible as the contact line is approached. Motivated by this observation, they dropped the viscous stress contribution in (1.5) leading to a “simplified GNBC”, formally reading as

$$-\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) = \sigma(\cos \theta_d - \cos \theta_e) \mathbf{n}_{\Gamma} \delta_{\Gamma}. \qquad (1.6)$$

It is evident that, by taking the inner product with the contact line normal vector, Eq. (1.6) can be formally reduced to an equation equivalent to (1.3) if the delta distribution is approximated with a regular function over a finite width. Indeed, Yamamoto et al. (2013, 2014) smoothed the delta distribution over a region of approximately four grid points. Using this estimate as the characteristic width of the delta function, the authors concluded that

$$U_{cl} \approx \frac{\lambda \sigma(\cos \theta_e - \cos \theta_d)}{\Delta \eta}, \qquad (1.7)$$

where $\Delta$ is the grid size. Obviously, the contact line speed in (1.7) can only be grid-independent if also the slip length is chosen in proportion to the grid size, i.e. if $\lambda \sim \Delta$. Consequently, they fixed the parameter $\chi := \lambda/\Delta$ in their simulations. The approach was extended by using the Cox-Voinov relation for $\theta_d$ in Yamamoto et al. (2014). Later, Yamamoto et al. (2016) used this method to study the withdrawing plate problem with a single wettable defect. Recently, the GNBC front-tracking approach was extended by Kawakami et al. (2023) using a so-called “rolling belt-model” inspired by the work of Lukyanov & Pryer (2017). Chen et al. (2019) used the GNBC in a Front Tracking method to study the coalescence-induced self-propelled motion of droplets on a solid surface. Shang et al. (2018) used a quite similar method to study droplet spreading and the motion of drops on surfaces subject to a shear flow. All these methods have in common that the uncompensated Young stress is distributed over a characteristic distance, which is related to the mesh size, see Ren & E (2007, 2011a,b); Zhang & Yue (2020).

In the present work, we propose a “sharp-interface, contact region GNBC” (CR-GNBC) formulation, where the contact line delta distribution is replaced by a smooth function with a characteristic width $\varepsilon > 0$. This width $\varepsilon$ is understood as a physical model parameter and is, therefore, chosen independently of the computational mesh. It has been shown recently by Kulkarni et al. (2023) that this model (i.e. the GNBC model with finite $\varepsilon$) admits a local $C^2$-regularity of the velocity in the vicinity of the moving contact line. We develop an implementation of the CR-GNBC in a geometrical Volume-of-Fluid method. This method turns out to be consistent with the fundamental kinematic law (1.4). For this purpose, the dynamic contact angle is not prescribed but is reconstructed from the volume fraction field in a neighborhood of the contact line. As one important preliminary step, we validate this “free contact angle” method by studying the advection problem by a prescribed velocity field (see Fricke et al. (2020)). In this case, the interface and the contact line are transported without a boundary condition for the contact angle and the results are validated against analytical solutions of (1.4). We couple this method to the CR-GNBC model and use the
<figure><img src="image_3.png" /><figcaption>Figure 1: Mathematical notation for the withdrawing tape setup.</figcaption></figure>

reconstructed contact angle $\theta_d$ as an input parameter to compute the uncompensated Young stress in the simulation.

## Structure of this article

We study the withdrawing tape problem as a prototypical example for a dynamic dewetting process. The setup follows the previous work by Afkhami et al. (2018). The solid wall is moving upwards with velocity $U_w \ge 0$ in the laboratory frame (see figure 1). So, we study the case of a receding contact line. We define the (global) capillary number with respect to the wall speed as

$$ Ca := \frac{\eta U_w}{\sigma}. \qquad (1.8) $$

For convenience, we define the contact line capillary number using the negative contact line speed, i.e. (note that the contact line speed $U_{cl}$ is always measured relative to the solid)

$$ Ca_{cl} := \frac{\eta(-U_{cl})}{\sigma}. \qquad (1.9) $$

In a quasi-stationary state, we have $-U_{cl} = U_w$ and hence $Ca = Ca_{cl}$. With this definition, we can always work with positive values for the capillary number. Notice that, in the literature, one will also find the convention that $Ca_{cl}$ is negative for a receding contact line and positive for an advancing contact line.

The mathematical derivation of the GNBC in a sharp-interface framework is revisited in Section 2. It is shown that the GNBC can be obtained as a combined closure relation for the dissipation due to slip along the liquid-solid surface and the contact line dissipation. Using the laws of kinematics, we derive the contact angle evolution equation in Section 2.4 and show that (1.3) holds for quasi-stationary states. Moreover, the GNBC thin film equation is derived in Section 2.5. The numerical implementation of the method in the geometrical Volume-of-Fluid solver is described in Section 3. We validate the numerical method by studying the kinematic transport of the contact angle and the curvature at the contact line. The results for the withdrawing tape problem are discussed in detail in Section 4. In particular, it is shown that the results are converging under mesh refinement. Notably, we
<table><caption>Table 1: List of Symbols</caption><thead><tr><th>Symbol</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>ρ</td><td>Density</td><td>kg/m<sup>3</sup></td></tr><tr><td>η</td><td>Viscosity</td><td>Pa·s</td></tr><tr><td>σ</td><td>Surface Tension</td><td>N/m</td></tr><tr><td><b>v</b></td><td>Velocity</td><td>m/s</td></tr><tr><td>β</td><td>Friction Coefficient</td><td>Pa·s/m</td></tr><tr><td><b>g</b></td><td>Gravitational Acceleration</td><td>m/s<sup>2</sup></td></tr><tr><td><b>D</b></td><td>Rate-of-Deformation Tensor</td><td>1/s</td></tr><tr><td><b>S</b> = 2η<b>D</b></td><td>Viscous Stress Tensor</td><td>Pa</td></tr><tr><td>Σ</td><td>Interface</td><td>-</td></tr><tr><td>V<sub>Σ</sub></td><td>Interface Normal Speed</td><td>m/s</td></tr><tr><td><b>n</b><sub>Σ</sub></td><td>Unit Interface Normal</td><td>-</td></tr><tr><td>κ = -∇<sub>Σ</sub> · <b>n</b><sub>Σ</sub></td><td>Interface Mean Curvature</td><td>1/m</td></tr><tr><td>Γ</td><td>Contact Line</td><td>-</td></tr><tr><td>U<sub>cl</sub></td><td>Contact Line Speed</td><td>m/s</td></tr><tr><td><b>n</b><sub>Γ</sub></td><td>Contact Line Normal (tangential to ∂Ω)</td><td>-</td></tr><tr><td>∂Ω</td><td>Solid Boundary</td><td>-</td></tr><tr><td><b>n</b><sub>∂Ω</sub></td><td>Unit Outer Normal to Ω</td><td>-</td></tr><tr><td>U<sub>w</sub></td><td>Wall Speed</td><td>m/s</td></tr><tr><td><b>U</b><sub>w</sub></td><td>Wall Velocity</td><td>m/s</td></tr><tr><td>λ = η/β</td><td>Slip Length</td><td>m</td></tr><tr><td>Ca</td><td>Wall Capillary Number</td><td>-</td></tr><tr><td>Ca<sub>cl</sub></td><td>Contact Line Capillary Number</td><td>-</td></tr><tr><td>Ca<sub>loc</sub></td><td>Capillary Number in the lab frame of reference</td><td>-</td></tr><tr><td>Ca<sub>tr</sub></td><td>Transition Capillary Number</td><td>-</td></tr><tr><td>ζ</td><td>Contact Line Friction</td><td>Pa·s</td></tr><tr><td>θ<sub>e</sub></td><td>Static Contact Angle</td><td>rad</td></tr><tr><td>θ<sub>d</sub></td><td>Dynamic Microscopic Contact Angle</td><td>rad</td></tr><tr><td>θ<sub>a</sub></td><td>Numerical angle one layer above the wall</td><td>rad</td></tr><tr><td>θ<sub>s</sub></td><td>Steady state contact angle</td><td>rad</td></tr><tr><td>t<sup>*</sup></td><td>Dimensionless time</td><td>-</td></tr><tr><td>t<sub>ε</sub></td><td>Slip length time scale</td><td>-</td></tr><tr><td>t<sub>r</sub></td><td>Relaxation time scale</td><td>-</td></tr></tbody></table>

can demonstrate by a mesh study that, unlike for the Navier slip model, the curvature at the contact line converges to a finite value. Finally, we conclude this study by an outlook to a non-linear variant of the GNBC, which can be derived as a non-linear closure based on the entropy production described earlier in Section 2.

## 2. Mathematical Modeling

### 2.1. Governing equations

We employ the sharp-interface continuum modeling approach. We start from the incompressible, two-phase Navier-Stokes equations with surface tension for Newtonian fluids under isothermal conditions (see, e.g., Slattery (1999); Prüss & Simonett (2016)). Inside the
fluid phases, the governing equations are

$$ \partial_t (\rho \mathbf{v}) + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) + \nabla p = \nabla \cdot \mathbf{S} + \rho \mathbf{g}, \quad (2.1) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad (2.2) $$

with the viscous stress tensor†

$$ \mathbf{S} = 2\eta\mathbf{D} = \eta(\nabla\mathbf{v} + \nabla\mathbf{v}^T). $$

These bulk equations are accompanied by jump conditions at the interface $\Sigma(t)$. The interface is modeled as a hypersurface (i.e. it has zero thickness) and separates the domain $\Omega$ into two bulk phases $\Omega_{1,2}(t)$ occupied by the two fluid phases (see figure 1). Assuming that no phase change occurs in the system, the normal component of the adjacent fluid velocities $\mathbf{v}_{1,2}$ at the interface are coinciding and equal to the speed of normal displacement $V_\Sigma$ of the interface, resulting in the kinematic boundary condition

$$ V_{\Sigma} = \mathbf{v} \cdot \mathbf{n}_{\Sigma} \quad \text{on} \quad \Sigma(t), \quad (2.3) $$

where $\mathbf{n}_\Sigma$ is the interface unit normal field. Additionally, no-slip between the fluid phases is usually assumed. Assuming further that the surface tension $\sigma$ is constant, the jump conditions for mass and momentum read as

$$ \llbracket \mathbf{v} \rrbracket = 0, \quad \llbracket p\mathbb{I} - \mathbf{S} \rrbracket \mathbf{n}_{\Sigma} = \sigma\kappa\mathbf{n}_{\Sigma} \quad \text{on} \quad \Sigma(t). \quad (2.4) $$

Here $\kappa := -\nabla_\Sigma \cdot \mathbf{n}_\Sigma$ is twice the mean curvature of the interface and

$$ \llbracket \psi \rrbracket (t, \mathbf{x}) := \lim_{h \to 0^+} (\psi(t, \mathbf{x} + h\mathbf{n}_\Sigma) - \psi(t, \mathbf{x} - h\mathbf{n}_\Sigma)) $$

is the jump of a quantity $\psi$ across the interface. We assume that the solid boundary $\partial\Omega$ is not able to store mass and we assume it to be impermeable. We consider an inertial frame of reference, where the wall is moving parallel to itself with a velocity $U_w \ge 0$ upwards (see figure 1). The impermeability condition in this frame of reference reads as

$$ \mathbf{v}_{\perp} = 0 \quad \text{on} \quad \partial\Omega, \quad (2.5) $$

where $\mathbf{v}_{\perp} = (\mathbf{v} \cdot \mathbf{n}_{\partial\Omega}) \mathbf{n}_{\partial\Omega}$ denotes the normal part of the velocity with respect to $\partial\Omega$.

In order to obtain a closed model, the system of equations (2.1)-(2.5) must be complemented by (one or more) additional boundary conditions describing

(i) the **wettability** of the solid (i.e. the static and dynamic contact angle) and

(ii) the **mobility** of the contact line (i.e. the tangential velocity $\mathbf{v}_\|$ at the solid boundary).

These boundary conditions are closure relations for the continuum mechanical description and must be thermodynamically consistent, i.e. they must obey in particular the second law of thermodynamics. To arrive at a consistent closure, we consider the available energy consisting of the kinetic energy of the bulk phases and the surface energies of the liquid-gas interface as well as the wetted area $W(t) \subset \partial\Omega$, i.e.

$$ E(t) := \int_{\Omega \setminus \Sigma(t)} \frac{\rho |\mathbf{v}|^2}{2} dV + \int_{\Sigma(t)} \sigma dA + \int_{W(t)} \sigma_w dA. $$

Here $\sigma = \sigma_{\text{lg}} > 0$ denotes the surface tension of the liquid-gas interface and

$$ \sigma_w = \sigma_{\text{ls}} - \sigma_{\text{sg}} $$

<footnote>† We use the symbol $\mathbf{D} = \frac{1}{2}(\nabla\mathbf{v} + \nabla\mathbf{v}^T)$ for the rate-of-deformation tensor.</footnote>
is the specific energy density for wetting the solid surface. Note that $\sigma_w$ is negative for a hydrophilic surface, as we see from Young's equation

$$ \sigma \cos \theta_e + \sigma_w = 0, \qquad (2.6) $$

which defines the “static” or “equilibrium” contact angle $\theta_e$. It is a purely mathematical exercise (see Ren & E (2007, 2011b), Fricke (2021) for details) to compute the rate of change $\dot{E}$ for a sufficiently regular solution of (2.1)-(2.5) (in the absence of external forces, i.e. for $\mathbf{g} = 0$). The result reads as

$$ \frac{dE}{dt} = - \int_{\Omega \setminus \Sigma(t)} \mathbf{S} : \mathbf{D} \, dV + \int_{\partial \Omega} (\mathbf{v}_{\|} - \mathbf{U}_w) \cdot (\mathbf{S} \mathbf{n}_{\partial \Omega})_{\|} \, dA + \int_{\Gamma(t)} \sigma (\cos \theta_d - \cos \theta_e) U_{\text{cl}} \, dl. \qquad (2.7) $$

In this formulation with a continuous velocity field, the scalar contact line speed (measured relative to the solid) is given as

$$ U_{\text{cl}} = \mathbf{v} \cdot \mathbf{n}_{\Gamma} - U_w. $$

Closure relations are required to satisfy the second law of thermodynamics† $\dot{E} \le 0$. The first contribution in (2.7) has a negative sign as we consider incompressible Newtonian fluids, i.e. $\mathbf{S} = 2\eta\mathbf{D}$. A linear closure for the second integral in (2.7) yields the well-known Navier slip condition, i.e.

$$ -\beta(\mathbf{v}_{\|} - \mathbf{U}_w) = (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\|} \quad \text{with a friction coefficient } \beta \ge 0. \qquad (2.8) $$

Using the slip length parameter $\lambda = \eta/\beta$, one may reformulate (2.8) as

$$ \mathbf{v}_{\|} + 2\lambda(\mathbf{D}\mathbf{n}_{\partial\Omega})_{\|} = \mathbf{U}_w. \qquad (2.9) $$

The third integral in (2.7) suggests that the dynamic contact angle $\theta_d$, which is mathematically defined as the angle of intersection‡ of the free surface $\Sigma$ with the solid boundary $\partial\Omega$, i.e.

$$ \cos \theta_d := -\mathbf{n}_{\Sigma} \cdot \mathbf{n}_{\partial\Omega} \quad \text{at } \Gamma(t), $$

should be linked to the contact line speed $U_{\text{cl}}$. A linear closure leads to the well-known condition

$$ -\zeta U_{\text{cl}} = \sigma(\cos \theta_d - \cos \theta_e) \quad \text{with } \zeta \ge 0, \qquad (2.10) $$

where $\zeta$ can be seen as a contact line friction coefficient. Note that also more general contact angle boundary conditions are possible if a non-linear closure relation is employed. To summarize, the “standard model”¶ based on the Navier slip condition is given by Equations (2.1)-(2.5), (2.9) together with (2.10) or a non-linear variant of the form

$$ \theta_d = f(U_{\text{cl}}) \quad \text{on } \Gamma(t). \qquad (2.11) $$

To ensure thermodynamic consistency, we require that

$$ \eta \ge 0, \quad \sigma \ge 0, \quad \lambda \ge 0, \quad U_{\text{cl}}(f(U_{\text{cl}}) - \theta_e) \ge 0. $$

† Note that we assume an isothermal system here. In this case, we may directly consider the change in available energy.

‡ Note that, in order to define the contact angle $\theta_d$, we have to assume that interface $\Sigma(t)$ has a well-defined normal field up to the boundary. This is the case even if the curvature has a logarithmic, hence integrable, singularity.

¶ The mathematical model (2.1)-(2.5), (2.9), (2.10) is one of the most commonly applied models for dynamic wetting in the literature. However, there are many more modeling approaches which aim at a regularization of the singularity and a prediction of the dynamics of wetting. For a survey of the field, we refer to the references de Gennes et al. (2004); Blake (2006); Shikhmurzaev (2008); Bonn et al. (2009); Snoeijer & Andreotti (2013a); Marengo & De Coninck (2022).
As shown by Fricke & Bothe (2020), for regular solutions, there is a fundamental inconsistency of boundary conditions in the standard model. In fact, the evolution of the contact angle is determined by the contact angle boundary condition (say (2.10)) as well as by the flow in the vicinity of the contact line according to (1.4). As a consequence, a regular solution of the system does not exist but a weak singularity is present at the contact line as shown already in Huh & Mason (1977).

## 2.2. Formal derivation of the GNBC

We now recall the derivation of the GNBC in the context of the sharp-interface framework described above starting from (2.7)

It is important to note that the GNBC was originally formulated in a *diffuse* interface framework (see Qian *et al.* (2003, 2006a)). However, the GNBC can be formally understood in the sharp interface model as a *combined closure* for the terms in the entropy production (2.7) that arise from the contact line motion and from slip at the solid-liquid boundary. The combined closure formally leads to a single boundary condition instead of the two independent conditions in the standard Navier slip model. Hence, the number of boundary conditions is reduced and one can show that the inconsistency at the contact line can be resolved this way (see below).

As a starting point, we consider the sum of the wall and the contact line dissipation, given as

$$ \mathcal{T} = \int_{\partial\Omega} (\mathbf{v}_{\parallel} - \mathbf{U}_w) \cdot (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} dA + \sigma \int_{\Gamma(t)} (\cos \theta_d - \cos \theta_e) U_{cl} dl. $$

By introducing the contact line delta distribution $\delta_{\Gamma}$, one can rewrite $\mathcal{T}$ as a single integral over the entire solid boundary $\partial\Omega$ according to

$$ \mathcal{T} = \int_{\partial\Omega} ((\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos \theta_d - \cos \theta_e) \mathbf{n}_{\Gamma}\delta_{\Gamma}) \cdot (\mathbf{v}_{\parallel} - \mathbf{U}_w) dA. $$

Note that it is possible to factor out the common co-factor $(\mathbf{v}_{\parallel} - \mathbf{U}_w)$ because the contact line speed can be written as $U_{cl} = (\mathbf{v}_{\parallel} - \mathbf{U}_w) \cdot \mathbf{n}_{\Gamma}$. A linear† closure relation is now provided by the generalized Navier boundary condition

$$ -\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) = (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos \theta_d - \cos \theta_e) \mathbf{n}_{\Gamma}\delta_{\Gamma} \quad \text{on } \partial\Omega \quad (2.12) $$

with a friction coefficient $\beta > 0$. Notice that the boundary condition given by (2.12), called “delta function GNBC” to distinguish it from other variants, should be understood in the sense of distributions.

There is an interesting link between (2.12) and the independent closure (2.9) together with (2.10): following Ren & E (2007), one may assume that the friction $\beta$ in (2.12) has a singular component, i.e.

$$ \beta = \beta_{\Gamma}\delta_{\Gamma} + \beta_{b} \quad \text{with} \quad \beta_{\Gamma}, \beta_{b} \ge 0. $$

In this case, equation (2.9) takes the form

$$ [\beta_{\Gamma}(\mathbf{v}_{\parallel} - \mathbf{U}_w) + \sigma(\cos \theta_d - \cos \theta_e) \mathbf{n}_{\Gamma}]\delta_{\Gamma} + [\beta_b(\mathbf{v}_{\parallel} - \mathbf{U}_w) + (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel}] = 0 \quad \text{on } \partial\Omega. $$

This distributional equation can be split into a singular and a regular part. This yields the

<footnote>† A non-linear generalization of the closure is discussed in Section 6.1.</footnote>
”standard closure relations” for slip and dynamic contact angle, i.e.

$$ \beta_b(\mathbf{v}_{\parallel} - \mathbf{U}_w) + (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} = 0, \quad \beta_{\Gamma}(\mathbf{v}_{\parallel} - \mathbf{U}_w) \cdot \mathbf{n}_{\Gamma} + \sigma(\cos \theta_d - \cos \theta_e) = 0. $$

Hence, the sharp-interface, sharp-contact line GNBC condition (2.12) and the standard closure (2.9) and (2.10) are formally equivalent if the friction $\beta$ has a singular component at the contact line. As we will see below, this singular component $\beta_{\Gamma}$ may arise as a ”lumped friction” that results from shrinking the contact zone to a line with vanishing thickness. Let us already note at this point, that we are not going to study the limit of a sharp contact line but model it as a finite contact region.

## 2.3. Contact Region GNBC (CR-GNBC) model

To obtain the CR-GNBC model, we replace the contact line delta distribution in (2.12) by a smooth function defined over a finite transition region with characteristic width $\epsilon$ such that

$$ \delta_{\Gamma}^{\epsilon} \ge 0, \quad \int \delta_{\Gamma}^{\epsilon}(x) \, dx = 1. $$

Note that this approach also requires to extend the definition of the contact angle $\theta_d$ and the contact line normal $\mathbf{n}_{\Gamma}$ away from the sharp contact line. In fact, the existence of the solid boundary, touching the interface at an angle strictly between 0 and $\pi$, provides a means for this extension of the contact line to a finite region. Then, the deviation of the contact angle from the equilibrium value appears in the velocity boundary condition leading to a force balance between sliding friction forces due to slip along the solid boundary, the tangential component of the viscous stress at the boundary and the uncompensated Young force:

$$ -(\mathbf{v}_{\parallel} - \mathbf{U}_w) = 2\lambda(\mathbf{D}\mathbf{n}_{\partial\Omega})_{\parallel} + \frac{\sigma}{\beta}[(\cos \theta_d - \cos \theta_e)\mathbf{n}_{\Gamma}\delta_{\Gamma}^{\epsilon}]. \quad (2.13) $$

Notably, the dynamic contact angle is not prescribed explicitly in this approach. Instead, the dynamics of the contact angle is determined by (2.13) and the kinematics of the interface transport.

## 2.4. Kinematics of the dynamic contact angle

We derive the evolution law for the contact angle, given a sufficiently regular solution of the CR-GNBC model (2.1)-(2.5) and (2.13). Below, we consider the limiting case of a free surface flow, where one phase is assumed to be a dynamically passive gas at a constant pressure. We assume that $\delta_{\Gamma}^{\epsilon}$ evaluated at the contact line yields a value of $1/\epsilon$. Then, the CR-GNBC condition, evaluated at the contact line, reads as

$$ \beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) + (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \frac{1}{\epsilon}\sigma(\cos \theta_d - \cos \theta_e)\mathbf{n}_{\Gamma} = 0 \quad \text{at } \Gamma. \quad (2.14) $$

By taking the inner product with the contact line normal vector $\mathbf{n}_{\Gamma}$ (normal to the contact line and tangential to the solid), we obtain the relation

$$ \frac{U_{cl}}{\lambda} + \langle \mathbf{n}_{\Gamma}, (\nabla \mathbf{v}) \mathbf{n}_{\partial\Omega} \rangle + \langle (\nabla \mathbf{v}) \mathbf{n}_{\Gamma}, \mathbf{n}_{\partial\Omega} \rangle + \frac{\sigma}{\epsilon\eta}(\cos \theta_d - \cos \theta_e) = 0 \quad \text{at } \Gamma. \quad (2.15) $$

Using the kinematic evolution equation for the contact angle derived in Fricke et al. (2019), one can show that the rate-of-change of the contact angle $\dot{\theta}_d$ is given as

$$ 2\dot{\theta}_d = - \langle \mathbf{n}_{\Gamma}, (\nabla \mathbf{v}) \mathbf{n}_{\partial\Omega} \rangle. \quad (2.16) $$

Moreover, it follows from the impermeability condition that the term $\langle \nabla \mathbf{v} \mathbf{n}_{\Gamma}, \mathbf{n}_{\partial\Omega} \rangle$ vanishes for a flat solid boundary. Therefore, we obtain the contact angle evolution law for a regular
solution of the CR-GNBC model. It reads as

$$ \dot{\theta}_d = \frac{U_{cl}}{2\lambda} + \frac{1}{\varepsilon} \frac{\sigma}{2\eta} (\cos \theta_d - \cos \theta_e). \qquad (2.17) $$

### Remarks

(i) Compared to the standard Navier slip model (see Fricke et al. (2019) for details), the uncompensated Young stress leads to an additional term in the equation for $\dot{\theta}_d$, which reads as

$$ \frac{1}{\varepsilon} \frac{\sigma}{2\eta} (\cos \theta_d - \cos \theta_e). $$

Obviously, the latter term is negative for $\theta_d > \theta_e$ (and positive for $\theta_d < \theta_e$) and, hence, drives the system towards equilibrium.

(ii) An important consequence of the CR-GNBC for quasi-stationary states is that it defines a functional dependence between the dynamic contact angle and the contact line speed. In fact, setting $\dot{\theta}_d = 0$ leads to the relation

$$ \mathrm{Ca}_{cl} = \frac{\eta(-U_{cl})}{\sigma} = \frac{\lambda}{\varepsilon} (\cos \theta_d - \cos \theta_e), \qquad (2.18) $$

or, equivalently, to

$$ -(\beta\varepsilon) U_{cl} = \sigma(\cos \theta_d - \cos \theta_e). \qquad (2.19) $$

By comparing (2.19) with (2.10), we see that the contact line friction parameter can be identified with the product of the “bulk friction” in the Navier slip condition and the width of the contact line region, i.e.

$$ \zeta = \beta\varepsilon. \qquad (2.20) $$

The latter equation has been proposed before by Blake et al. (2015) in the context of the molecular kinetic theory. Physically, it indicates that, within the present modeling framework, there is only one friction mechanism that affects both the slip at the solid boundary and the dynamics of the microscopic contact angle. In this sense, the contact line friction $\zeta$ can be understood as the lumped wall friction of the contact region.

(iii) From (2.16) we conclude that the stress component $\langle \mathbf{n}_\Gamma, (\nabla \mathbf{v}) \mathbf{n}_{\partial\Omega} \rangle$ vanishes at the contact line for quasi-stationary states, i.e. for $\dot{\theta} = 0$. Hence, there appears to be “perfect slip” at the contact line in that case. Actually, the concepts of the “apparent slip length” $\lambda_a$ (see figure 2) and the physical slip parameter defined as $\lambda = \eta/\beta$ must be distinguished for the GNBC model. In fact, the uncompensated Young stress is able to reverse the sign of the velocity gradient at the contact line. In this case, fluid particles at the solid boundary may have a larger tangential velocity than fluid particles slightly above the boundary. This situation corresponds to a negative apparent slip length (see figure 2). It is, however, caused by the uncompensated Young stress in the velocity boundary condition. The physical slip length parameter $\lambda$ is still positive and finite in all cases.

(iv) Note that (2.17) can be rephrased as a generalized mobility law of the form

$$ U_{cl} = f(\theta_d, \dot{\theta}_d). \qquad (2.21) $$

Therefore, the contact line speed depends on the contact angle $\theta_d$ but also on its rate-of-change $\dot{\theta}_d$ which, in turn, can be computed from $\nabla \mathbf{v}$ (see Fricke et al. (2019)). In this sense, the contact line speed in the GNBC model depends on the flow field in the vicinity of the contact line.

(v) Moreover, the GNBC can be understood as an inhomogeneous Robin condition for
<figure><img src="image_4.png" /><figcaption>Figure 2: Different cases for the apparent slip length $\lambda_a$: positive, perfect and negative slip (reference frame with $U_w = 0$).</figcaption></figure>

the velocity. Hence, the GNBC enforces a flow whenever $\theta_d \neq \theta_e$. In contrast to the standard Navier slip model, the GNBC model is able to describe the relaxation process of the contact angle.

## 2.5. The CR-GNBC thin film equation

We now derive the thin film equation for the CR-GNBC and compare it with other known thin film equations in the context of dynamic contact lines. The coordinate system is shown in figure 1. Under the thin film assumption, we consider that the pressure remains constant along the y-axis and that the Laplace pressure jump across the interface can be expressed as

$$ \Delta p = - \frac{\sigma H''}{(1 + H^2)^{3/2}} \approx -\sigma H''. $$

The x-momentum equation

$$ -\frac{1}{\rho} \frac{\partial p}{\partial x} + g + \frac{\eta}{\rho} \frac{\partial^2 v}{\partial y^2} = 0 \qquad (2.22) $$

is supplemented by a free surface condition

$$ \frac{\partial v}{\partial y}(x, H) = 0. $$

Since we are in the reference frame of the contact line, we can write the CR-GNBC as

$$ v(x, 0) + \lambda \frac{\partial v}{\partial y}(x, 0) = U_w \left[ 1 - f\left(\frac{x}{\epsilon}\right) \right], \qquad (2.23) $$

where $f\left(\frac{x}{\epsilon}\right)$ is the smoothed Dirac function that will be defined in (3.6). For further details on the formulation in the reference frame of the contact line, we refer to Kulkarni et al. (2023). Given the contact line boundary condition, the free surface condition and the Laplace pressure jump, the x-momentum equation (2.22) can now be written as an ordinary differential equation in terms of $H$

$$ H''' + \frac{1}{l_c^2} = - \frac{3 \mathrm{Ca} \left(1 - f\left(\frac{x}{\epsilon}\right)\right)}{H(H - 3\lambda)} - \frac{-3\eta Q}{\sigma H^2 (H - 3\lambda)}, \qquad (2.24) $$

where $l_c$ is the capillary length and $Q = - \int_0^H v \, ds$ is the total flux. Assuming steady state, where $Q = 0$ and $H''' \gg 1/l_c^2$ in the vicinity of the contact line, we obtain the CR-GNBC thin film equation

$$ H''' = \frac{3 \mathrm{Ca} \left( \tanh^2 \left( \frac{x}{\epsilon} \right) \right)}{H(H - 3\lambda)}. \qquad (2.25) $$
<table><caption>Table 2: Thin film equations for various contact line boundary conditions.</caption><thead><tr><th>Boundary condition</th><th>Thin film equation form</th><th>Smoothness</th></tr></thead><tbody><tr><td>No slip: v = U<sub>w</sub></td><td>H''' = <sup>1</sup>&frasl;<sub>H<sup>2</sup></sub><br>Duffy &amp; Wilson (1997)</td><td>Singular H'<br>(angle singularity)</td></tr><tr><td>Navier slip: v &minus; &lambda;<sup>∂v</sup>&frasl;<sub>∂y</sub> = U<sub>w</sub></td><td>H''' = <sup>3Ca</sup>&frasl;<sub>H<sup>2</sup> + 3&lambda;H</sub><br>Eggers (2004)</td><td>Singular H'' &sim; log x<br>(curvature singularity)</td></tr><tr><td>Super-slip: v &minus; &lambda;<sup>2</sup><sup>∂<sup>2</sup>v</sup>&frasl;<sub>∂y<sup>2</sup></sub> = U<sub>w</sub></td><td>H''' = <sup>Ca</sup>&frasl;<sub>H<sup>2</sup> + &lambda;<sup>2</sup></sub><br>Hocking (2001)</td><td>Regular H''<br>(singularity-free)</td></tr><tr><td>Super-slip: v &minus; &lambda;<sup>∂v</sup>&frasl;<sub>∂y</sub> &minus; &lambda;<sup>2</sup><sup>∂<sup>2</sup>v</sup>&frasl;<sub>∂y<sup>2</sup></sub> = U<sub>w</sub></td><td>H''' = <sup>Ca</sup>&frasl;<sub>H<sup>2</sup>/3 + &lambda;H + &lambda;<sup>2</sup></sub><br>Devauchelle et al. (2007)</td><td>Regular H''<br>(singularity-free)</td></tr><tr><td>CR-GNBC</td><td>H''' = <sup>3Ca(tanh<sup>2</sup><sup>x</sup>&frasl;<sub>&epsilon;</sub>)</sup>&frasl;<sub>H(H &minus; 3&lambda;)</sub><br>current paper</td><td>Regular H''<br>(singularity-free)</td></tr></tbody></table>

From (2.25), we can see that for $x \ll \epsilon$, $H''$ does not diverge and approaches a constant value at the contact line ($x = 0$). Hence, the equation is singularity-free. Our CR-GNBC model can therefore be viewed as Navier slip with a *smoothening well* of width $\epsilon$ around the contact line where the uncompensated Young stress acts. A comparison of thin film equations from the literature and their respective smoothness is presented in Table 2.

### 3. Numerical Methods

#### 3.1. *The Volume-of-Fluid method*

The Volume-of-Fluid (VOF) method for representing fluid interfaces coupled with a flow solver is well-known to be suited for solving interfacial flows (see e.g. Scardovelli & Zaleski (1999); Popinet & Zaleski (1999); Tryggvason et al. (2011); Marić et al. (2020)). We use the free software *Basilisk*, a platform for the solution of partial differential equations on adaptive Cartesian meshes (Popinet (2009, 2015, 2018)). For a two-phase flow, the volume fraction $c(\mathbf{x}, t)$ is defined as the integral of the first fluid's characteristic function in the control volume. The volume fraction $c(\mathbf{x}, t)$ is used to define the density and viscosity in the control volume according to

$$ 
\begin{aligned}
\rho(c) &\equiv c\rho_1 + (1-c)\rho_2, \\
\mu(c) &\equiv c\mu_1 + (1-c)\mu_2,
\end{aligned} 
\qquad (3.1) 
$$

with $\rho_1$, $\rho_2$ and $\mu_1$, $\mu_2$ the densities and viscosities of the phase 1 and 2 respectively. The advection equation for the density is then replaced by the equation for the volume fraction

$$ \partial_t c + \mathbf{v} \cdot \nabla c = 0. \qquad (3.2) $$
The projection method is used to solve the incompressible Navier-Stokes equations combined with a Bell-Collela-Glaz advection scheme and a VOF method for interface tracking. The resolution of the surface tension term is directly dependent on the accuracy of the curvature calculation. The height-functions method, described in Afkhami & Bussmann (2008, 2009), 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 grid direction closest to the normal of the interface. In two dimensions, a 7 × 3 stencil around an interface cell is constructed and the heights are evaluated by summing volume fractions horizontally, i.e.

$$h_j = \sum_{k=i-3}^{k=i+3} c_{j,k} \Delta, \qquad (3.3)$$

with $c_{j,k}$ the volume fraction and $\Delta$ the grid spacing. The heights are then used to compute the interface normal $\mathbf{n}_{\Sigma}$ and the curvature $\kappa$ according to

$$\mathbf{n}_{\Sigma} = (h_x, -1),$$

$$\kappa = \frac{h_{xx}}{(1 + h_x^2)^{3/2}}, \qquad (3.4)$$

where $h_x$ and $h_{xx}$ are discretized using second-order central differences. The orientation of the interface, characterized by the contact angle – the angle between the normal to the interface at the contact line and the normal to the solid boundary – is imposed in the contact line cell. It is important to note that a numerical specification of the contact angle affects the overall flow calculation in two ways:

(i) it defines the orientation of the interface reconstruction in cells that contain the contact line;

(ii) it influences the calculation of the surface tension term by affecting the curvature computed in cells at and near the contact line.

We now present the numerical implementation of the Generalized Navier Boundary Condition as written in (2.14). The boundary condition is applied on the solid surface with a smoothing function that takes into account the relative position along the boundary with respect to the contact line, denoted $x$

$$\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) + (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + f\left(\frac{x}{\varepsilon}\right) \sigma(\cos\theta_d - \cos\theta_e) \mathbf{n}_{\Gamma} = 0 \quad \text{on } \partial\Omega, \qquad (3.5)$$

with $f\left(\frac{x}{\varepsilon}\right)$ the smoothed Dirac function defined as

$$f\left(\frac{x}{\varepsilon}\right) = \frac{\left(1 - \tanh^2\left(\frac{x}{\varepsilon}\right)\right)}{\varepsilon}. \qquad (3.6)$$

This specific function is smooth, symmetric and preserves the area for varying $\varepsilon$, characteristic features that are necessary for the well-posedness of the discrete boundary condition. The boundary condition can be expressed as an inhomogeneous Robin boundary condition for the parallel velocity $\mathbf{v}_{\parallel}$, as outlined above:

$$\mathbf{v}_{\parallel} + \frac{1}{\beta}(\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} = \mathbf{U}_w + \frac{1}{\beta}f\left(\frac{x}{\varepsilon}\right) \sigma(\cos\theta_e - \cos\theta_d) \mathbf{n}_{\Gamma} \quad \text{on } \partial\Omega. \qquad (3.7)$$

We use the Navier boundary condition (Navier slip) that was implemented in the same framework in Fullana et al. (2020) and tested as a localized slip boundary condition in Lācis
<figure><img src="image_5.png" /><figcaption>Figure 3: Extrapolation of the contact angle $\theta_d$ using the angle above $\theta_a$ located $3/2 \Delta$ away from the wall. $h_0$ to $h_2$ denote the horizontal heights.</figcaption></figure>

et al. (2020). The difference lies now in the space dependent right-hand side of (3.7). The uncompensated Young’s stress that only acts at the contact line through the discrete Dirac function, needs to be computed at each grid point.

The numerical approach in this study stands out for its free contact angle method. Instead of setting the dynamic angle $\theta_d$, we reconstruct it from the interface geometry and use it as an input parameter to calculate the right-hand side of (3.7). To reconstruct such a consistent angle from the volume fraction field, we use a Taylor expansion of the contact angle along the coordinate direction normal to the boundary (see y-axis in figure 3)

$$ \theta_d = \theta_a + \delta \frac{d\theta}{dy} + O(\delta^2), \qquad (3.8) $$

where $\theta_a$ is the angle above located at the distance $y = \delta$ from the wall. Using that $d\theta/ds = \kappa$ and projecting onto the wall-normal direction, we obtain

$$ \theta_d = \theta_a + \delta \frac{\kappa \sqrt{1+x'}}{\sin \theta_a} + O(\delta^2), \qquad (3.9) $$

with $x'$ is the local slope of the interface. Discretely, by setting $\delta = 3/2 \Delta$, such that the angle above is computed one layer above the wall and using height-functions to compute the slope of the interface, such that $h_y$ is the first order derivative of the height-functions in the y direction (normal to the wall), we obtain

$$ \theta_d = \theta_a + \frac{3}{2} \Delta \frac{\kappa \sqrt{1 + h_y^2}}{\sin \theta_a} + O(\Delta^2). \qquad (3.10) $$

Figure 3 provides a schematic illustration of this extrapolation process. Once the extrapolated angle is computed, we enforce it through appropriate local modification of height-functions in the ghost layer, similar to a regular contact angle. Algorithm 1 is a concise summary of the two-step procedure to apply the CR-GNBC in the VOF framework.

## 3.2. Kinematic transport of the contact angle

We validate the free contact angle method presented in (3.10) through an analysis of the kinematic transport of the contact angle in a simplified setup. Leveraging kinematic considerations, Fricke et al. (2020) and Fricke (2021) derived analytical solutions for the
<caption>**Algorithm 1:** CR-GNBC pseudo-code</caption>

**for each boundary cell do**
1. Locate the contact line cell
2. Locate the cell one grid point above the contact line
3. Compute the angle above $\theta_a$ using the unit normal $\mathbf{n}_\Sigma$
4. Compute the first order derivative $h_x$ of the height-functions
5. Compute the interface curvature $\kappa$
6. Compute the extrapolated angle $\theta_d$ using (3.10)
**end**
7. Apply $\theta_d$ at the contact line through height-functions

**for each boundary cell do**
8. Compute the right-hand-side of (3.7) using $\theta_d$ and $\theta_e$
**end**
9. Apply the boundary condition for $\mathbf{v}_\parallel$ using (3.7)

transport of the contact angle and the curvature for some specific velocity fields. To validate the present approach within the VOF framework, we conduct advection testcases for an initially circular interface in contact with the domain boundary. These advection tests are carried out for various grid sizes.

The setup involves a disk with a dimensionless diameter $D = 1$ in a $2 \times 2$ domain. The center of the disk is shifted by $\delta_s = 0.05$ above the substrate, resulting in an initial contact angle of $\theta_0 = \arccos(\delta_s/R) \approx 84.26^\circ$. The velocity field across the entire domain is defined as:

$$ 
\begin{aligned}
v_x &= c_1 \cos(\pi t) x + c_2 \cos(\pi t) y, \\
v_y &= -c_1 \cos(\pi t) y.
\end{aligned} 
\qquad (3.11) 
$$

Here, $v_x$ and $v_y$ represent the x and y components of the velocity, while $c_1$ and $c_2$ are positive constants. We aim to validate the accuracy and reliability of the angle extrapolation method under varying grid sizes, where we only consider the advection equation of the color function (3.2).

The prescribed incompressible velocity field (3.11) will induce oscillations of the interface in both vertical and horizontal directions. The angle formed at the contact line is determined by this motion and varies in time. From the relations derived in Fricke et al. (2020), we compare the observed numerical contact angle with the analytical value $\theta_{an}$, given by the formula

$$ \theta_{an}(t) = \frac{\pi}{2} + \tan^{-1} \left( \frac{-1}{\tan \theta_0} e^{2c_1 S(t)} + \frac{c_2}{2c_1} (e^{2c_2 S(t)} - 1) \right) \qquad (3.12) $$

with

$$ S(t) = \frac{\sin(\pi t)}{\pi}. \qquad (3.13) $$

Moreover, we validate the evolution of the curvature by comparison with the reference one $\kappa_{an}$, which is given as the solution of the ordinary differential equation (see Fricke (2021))

$$ \frac{d \kappa_{an}}{dt} = -3 \kappa_{an} \cos(\pi t) [c_1 \cos^2(\theta_{an}) - c_2 \cos(\theta_{an}) \sin(\theta_{an}) - c_1 \sin^2(\theta_{an})] \qquad (3.14) $$

with the initial condition $\kappa_0 = 2/D = 2$. We use Vofi (see Bnà et al. (2016)) to initialize the volume fraction field, which yields a more accurate initial curvature. We conduct simulations
<figure><img src="image_6.png" /><figcaption>Figure 4: Validation of the free angle extrapolation method for varying grid sizes. (a) Temporal evolution of the contact angle. (b) Temporal evolution of the curvature. In both panels, blue corresponds to $D/\Delta = 32$, orange to $D/\Delta = 64$, green to $D/\Delta = 128$, red to $D/\Delta = 256$, and the black dashed line is the analytical solution. (c) Convergence of both the angle and curvature errors with $L_2$ and $L_\infty$ norms, showing quasi-second-order convergence for the angle and quasi-first-order convergence for the curvature.</figcaption></figure>

with $c_1 = 0.5$ and $c_2 = 0.2$, using the free angle extrapolation method. The simulations run until a final dimensionless time $T = 10$, and we examine the convergence of the method with grid sizes varying from 32 to 256 points per diameter (see figure 4) with a fixed timestep $\delta t = 0.2 \Delta/\max(|c_1|, |c_2|)$. The extracted contact angles and curvatures show convergence towards the analytical solution. As shown in figure 4(c), we observe quasi-second-order convergence for the contact angle and quasi-first-order convergence for the curvature, the latter is expected since $\kappa \sim \partial\theta/\partial s$, with $s$ the arc-length. We emphasize that this validation corresponds to a purely kinematic advection problem, in which convergence is expected to degrade by one order from interface position to contact angle and to curvature, since each quantity involves an additional spatial derivative of the interface geometry. Such behavior is specific to purely advective test cases, while in fully coupled dynamic simulations, surface-tension-induced regularization typically improves curvature convergence (Popinet 2009, 2018).
<figure><img src="image_7.png" /></figure>

Figure 5: Steady-state meniscus example for Ca = 0.1 using the present CR-GNBC with $\epsilon = 0.05$. The image is in the contact line's reference frame, where the left plate is pulled up with $\bar{U}_w = \sqrt{Ca}$. The inset shows a zoom around the contact line with streamlines highlighting a stagnation point in the upper phase

# 4. Results

We apply the numerical method for the CR-GNBC model to the pulling plate setup, following the approach discussed in Section 1. This setup is akin to setup B of Afkhami et al. (2018). We set the density ratio to $\rho_1/\rho_2 = 5$ and the viscosity ratio to $\mu_1/\mu_2 = 1$. The dimensionless capillary length is set to $l_c = 1$, the dimensionless viscosities to $\bar{\mu}_1 = \bar{\mu}_2 = \bar{U}_w$ and the dimensionless density $\bar{\rho}_1 = 1$, such that the Reynolds number is unity ($\text{Re} = \bar{\rho}_1 \bar{l}_c \bar{U}_w / \bar{\mu}_1 = 1$) for any prescribed wall capillary number. In this configuration, setting the dimensionless surface tension to $\bar{\sigma} = 1$, the wall velocity is $\bar{U}_w = \sqrt{Ca}$.

Figure 5 displays the results of a steady-state simulation. The image is presented in the reference frame of the contact line, where the contact line remains stationary while the left wall is pulled upwards. The velocity field relaxes, creating a stagnation point at the contact line. Additionally, the streamlines reveal another stagnation point formed above the interface in the lighter phase. It is worth noting that the characteristics of this additional stagnation point depend on the viscosity ratio, although our primary focus is not on this aspect.

In the pulling plate setup, a distinctive characteristic is the presence of a de-wetting transition capillary number $Ca_{tr}$, marking the point beyond which liquid film entrainment occurs, leading to an absence of a steady-state position for the contact line. Previous numerical results by Afkhami et al. (2018) identified this transition capillary number, but it was grid-dependent. Using the CR-GNBC method, with $\epsilon$ resolved (i.e. larger than the grid size $\Delta$), we obtain a grid-independent $Ca_{tr}$. This is depicted in figure 6, which shows the contact line position representing the fluid film height over time. For $Ca \le 0.12$, a steady-state height is eventually reached; however, for $Ca = 0.14$, the height continually increases. Thus, we determine that $Ca_{tr}$ for this case is $Ca = 0.13 \pm 0.01$. The influence of Young's stress is evident when comparing figure 6(a) with figure 6(b). The $Ca_{tr}$ remains the same, but the steady-state height exhibits a slight decrease. A convergence study demonstrating the grid independence of the CR-GNBC is presented in Appendix A. Furthermore, with the present
<figure><img src="image_8.png" /><figcaption>Figure 6: Vertical height of the contact line as a function of time for different capillary numbers Ca, presented separately for (a) simple Navier boundary condition and (b) CR-GNBC. Steady-state heights are achieved, and a transition $Ca_{tr}$ is observed, beyond which the liquid film rises continuously. In both (a) and (b), $Ca_{tr} = 0.13$. Simulations are conducted with $\epsilon = 0.05$, $\theta_c = 90^\circ$, and a resolution of $\epsilon/\Delta = 5.12$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: (a) Evolution of the dynamic contact angle $\theta_d$ in the CR-GNBC simulation for various Ca. The angle begins to deviate from the initial value of 90° and eventually reaches a steady state. Around $Ca_{tr}$, the angle exhibits oscillations over time. (b) The relaxation plot on a $\theta$ – Ca plane. Here, $Ca_{loc}$ represents the contact line capillary number. Time progresses from right to left, and a maximum in $Ca_{loc}$ is reached at $t_\epsilon = 1$, which corresponds to the slip length timescale ($\epsilon/U_w$). After this point, $Ca_{loc}$ starts relaxing towards a steady state ($Ca_{loc} = 0$). Above $Ca_{tr}$, $Ca_{loc}$ reaches a minimum and starts rising again. This set of simulations are the same as in figure 6(b).</figcaption></figure>

approach, the parameters influencing $Ca_{tr}$ are $\epsilon$, the slip length $\lambda$ and the equilibrium contact angle $\theta_c$. The dependence of these parameters on the $Ca_{tr}$ is presented in Appendix B.
4.1. Relaxation towards steady state

Starting from a horizontal two-fluid interface at rest, we now compare the transient characteristics. A distinctive feature of the present CR-GNBC model is that the contact angle is not fixed a priori. Figure 7(a) illustrates the contact angle $\theta_d$ as a function of time. The angle initiates at 90° and subsequently relaxes to a steady-state value different from 90°. Despite converging to a steady state, the observed value of $\theta_d$ exhibits spurious oscillations. These oscillations intensify with increasing Ca; however, their influence is minor, with amplitudes remaining below 0.5° and diminishing with grid refinement. In figure 7(a), we observe an interesting trend when plotting $\theta_d$ against $Ca_{loc}$, as shown in figure 7(b). Here, $Ca_{loc}$ represents the contact line $Ca$ in the lab frame. It starts at 0 since everything is initially at rest and eventually returns to 0 in a quasi-stationary state. During the transient phase, although we set the solid velocity to $U_w$ instantly, $Ca_{loc}$ takes some time to reach its maximum value. This time, defined as $t_\epsilon = \epsilon/U_w$, represents a relaxation timescale due to contact line friction. While a detailed examination of the behavior for $t < t_\epsilon$ is beyond the current study's scope, we observe that, once $Ca_{loc}$ reaches its peak, it begins to relax to the steady state where $Ca_{loc} = 0$, and $\theta_d$ follows the GNBC law (2.18).

In figure 8, we illustrate the behavior of each term of the CR-GNBC equation (3.7). We analyze and present each outcome for three different boundary conditions:

(i) Navier slip with a constant contact angle $\theta_d = \theta_e$:

$$v_{\parallel} + \frac{1}{\beta}(Sn\partial\Omega)_{\parallel} = U_w \quad \text{on } \partial\Omega, \qquad (4.1)$$

(ii) No slip with uncompensated Young stress, with the "free angle" method (3.10):

$$v_{\parallel} = U_w + \frac{1}{\beta}f\left(\frac{x}{\epsilon}\right)\sigma(\cos\theta_e - \cos\theta_d) n_{\Gamma} \quad \text{on } \partial\Omega, \qquad (4.2)$$

(iii) Full CR-GNBC as written in (3.7) which combines contributions from both the above cases.

We conducted simulations for each individual case (i), (ii), and (iii) and illustrate the behavior of each term in figure 8. In figure 8(a) we see the angle as a function of time. Since we start from a horizontal surface, all plots begin at 90°. The green curves, representing the Navier slip case, converge to the constant imposed value of 90°. The relaxation to a steady-state angle is accompanied by oscillations, whose amplitude decreases with grid refinement. The blue curves, representing the behavior of uncompensated Young's stress with a no-slip boundary condition, show the effect of the free contact angle. Because the Young stress term involves the free contact angle method, the steady-state angle differs from 90° and relaxes to the GNBC law contact angle as the grid is refined. These spurious oscillations are less pronounced than in the Navier slip case. Finally, the red curves represent the full CR-GNBC model. At the same level of grid refinement, the CR-GNBC model outperforms the Young stress case (blue curves) by being closer to the expected GNBC law contact angle and outperforms the Navier slip (green curves) by having fewer spurious oscillations. In the plot of $Ca_{loc}$ vs time, we see that, although all the curves eventually relax to the steady state of $Ca_{loc} = 0$, there is a difference in the initial relaxing stage. As soon as the simulation is started we see that since blue curves have no slip, they rise to the $Ca_{loc} = Ca$ in dimensionless time $t_\epsilon$ and then relax to the steady state value, while the CR-GNBC and slip cases rise to the value equal to $Ca_{loc} < Ca$.

Figure 8(b) shows the phase diagram on a $Ca_{loc}-\theta$ plane. This figure sums up the overall behaviour of the contact line dynamics in each case and a characteristic behaviour of each set could now be identified. The timeline in this figure progresses from right to left.

(i) In the Navier slip case, we observe that at $t = 0$ and for $\theta_d = 90°$, when the interface
<figure><img src="image_10.png" /><figcaption>Figure 8: The relaxation plots Navier slip (green curves), CR-GNBC (red curves) and no-slip with Young stress (blue curves). All the plots are done for Ca = 0.04 and ε = 0.05. The grid resolution is reported in terms of ε/Δ, and color intensity is increased to show higher resolution. Figure (a) shows the contact angle θd and the contact line speed in the lab frame of reference Ca<sub>loc</sub> as a function of time. In (b) we show the phase diagram resulting from figure (a). The dashed black line represents the GNBC law angle in the steady state. Time flows from right to left and aligns the curves. Each curve set has its own characteristic feature. The oscillations, present in (a), are faded in the phase diagram for clarity.</figcaption></figure>

is horizontal, Ca<sub>loc</sub> is null. Then, Ca<sub>loc</sub> suddenly rises to a maximum value, which remains lower than the imposed Ca. This rapid rise occurs within the relaxation time t<sub>ε</sub>, where ε is the slip length. This behavior aligns with the discussion in figure 7(b). Subsequently, the contact line relaxes to a steady state where Ca<sub>loc</sub> returns to zero. This relaxation is accompanied by spurious oscillations in the contact angle θ<sub>d</sub>. Ideally, in this case, the system should relax to θ<sub>d</sub> = 90° throughout the motion and also in the steady state (given that we impose a constant θ<sub>d</sub> = θ<sub>e</sub> = 90°), which is indeed observed as the grid is refined. The final angle θ<sub>d</sub> converges to 90°, and spurious oscillations diminish with increasing grid refinement.
<figure><img src="image_11.png" /><figcaption>Figure 9: Curvature profiles relative to the radial distance from the contact line. The red curves represent curvature under the Navier slip boundary condition (slip), showing a logarithmic divergence. In contrast, the blue curves (GNBC) demonstrate the convergence to a finite curvature value and, thus, the removal of the singularity present in the NBC. Simulations are conducted with Ca = 0.08 and ε = 0.05. The equilibrium angle is θₑ = 90°, and Δ denotes the grid size. Various color intensities denote grid refinement, where lighter shades correspond to a coarse mesh, and darker shades indicate a fine mesh.</figcaption></figure>

(ii) In the no-slip with Young's stress, we notice an interesting pattern. At the start (t=0), the simulation begins with Ca<sub>loc</sub> = 0 and θ<sub>d</sub> = 90° at the lower right of figure 8. However, as soon as we advance in time, Ca<sub>loc</sub> increases to a maximum value equal to Ca, and then, starts relaxing to 0. With the presence of uncompensated Young stress, it ideally should relax to the GNBC law contact angle indicated by the dashed line in figure 8. We observe that oscillations are decreasing with grid refinement and the final value of the contact angle is converging towards the GNBC law angle.

(iii) In the CR-GNBC case, we observe characteristics from both (i) and (ii). Initially, both Ca<sub>loc</sub> and θ<sub>d</sub> start from zero. Subsequently, Ca<sub>loc</sub> reaches a maximum during the relaxation time and eventually relaxes to the GNBC law contact angle. The notable advantage of the CR-GNBC is that even with a modest resolution of 5 grid points per slip length, the spurious oscillations, compared to case (i) at the same resolution, are significantly reduced. Moreover, the accuracy in relaxing towards the GNBC law contact angle (dashed line) is substantially improved compared to case (ii). Further grid refinement leads to a continued reduction in spurious oscillations and enhances accuracy.

## 4.2. Steady-state contact line dynamics: the GNBC smoothing signature

We now demonstrate the full regularization of the contact line singularity achieved by the present CR-GNBC method. Figure 9 presents the curvature as a function of the distance from the contact line for various grid resolutions. The Navier slip model exhibits a logarithmic divergence in curvature, consistent with the analytical findings of Devauchelle et al. (2007) and Kulkarni et al. (2023). While the singularity in the Navier slip model is integrable and considered “weak”, it induces a pressure singularity, rendering the slip model physically ill-posed. In contrast, the present CR-GNBC model regularizes the logarithmically singular curvature at the contact line (κ ~ log r), establishing it as a physically well-posed model.

In Section 2, we showed that assuming a C<sup>1</sup> velocity field up to the contact line in the reference frame of the moving wall, the rate of change of the contact angle scales with the
<figure><img src="image_12.png" /><figcaption>Figure 10: (a) Wall shear stress in a steady-state simulation plotted against the vertical position y. The dashed line corresponds to the Navier slip boundary condition, while the solid line corresponds to the CR-GNBC case. Both curves largely overlap, except for a small region shown in the zoom-ins for Navier slip and CR-GNBC in (b). The zoomed-in figures are normalized by the contact line position, where 0 on the x-axis corresponds to the contact line position. Notably, in the CR-GNBC case, the shear stress at the contact line is zero, whereas it is not the case for the Navier slip. The simulations are conducted with fixed Ca = 0.08 and ε = 0.05 and varying grid sizes.</figcaption></figure>

shear stress at the contact line. In steady state, where $\dot{\theta}_d = 0$, the shear stress must approach zero as it reaches the contact line. A non-zero shear stress would indicate a violation of the smoothness assumption made by Fricke et al. (2019). This violation occurs in the Navier slip model, as non-zero shear stress is necessary for contact line motion. In figure 10, we observe the behavior of shear stress for both Navier slip and CR-GNBC in steady state. At the contact line, the shear stress converges to zero within the $\epsilon$ region in the CR-GNBC case, aligning with the expected smoothness of the flow field. However, for the Navier slip model, the shear stress fails to converge to zero.

Having demonstrated that the shear stress at the contact line in the steady state is zero using
<figure><img src="image_13.png" /><figcaption>Figure 11: The behavior of the quasi-stationary value of $\theta_d$ vs. Ca is illustrated for various $\theta_e$ and compared with the GNBC law (2.18). The solid lines represent the analytical expression of the steady-state behavior expected from (2.18), while the dots depict simulation results. The different colors represent various $\theta_e$, progressing from left to right (black to red) as 45°, 60°, 75°, 90°, and 120°, respectively. The horizontal lines denote the $Ca_{tr}$ for each equilibrium angle considered. An excellent agreement between simulations and the GNBC law is observed up to $Ca < Ca_{tr}$.</figcaption></figure>

the CR-GNBC, we proceed to compare the quasi-stationary state GNBC relation (2.18) with our simulation results in figure 11. Remarkably, we observe excellent agreement between the simulation outcomes and the quasi-stationary GNBC law, particularly for $Ca < Ca_{tr}$. It is essential to note that the behavior of $Ca_{cl} = f(\theta_s)$ in figure 11, as predicted by the quasi-stationary GNBC law (2.15), is not explicitly imposed but is a direct outcome from the simulations.

# 5. Conclusion

To summarize, we have developed an implementation of the Contact Region Generalized Navier Boundary Condition (CR-GNBC) in a geometrical Volume-of-Fluid method. In this method, the dynamic contact angle is not prescribed but is controlled by kinematics through the velocity boundary condition. This is achieved by reconstructing the contact angle at the boundary using the interface normal and the curvature one cell layer away from the boundary. We validate the resulting free contact angle method by studying the interface advection problem in the presence of a moving contact line in Section 3.2. In the present approach, the uncompensated Young stress is distributed over a characteristic width $\epsilon$, that is defined independently of the mesh size. Using the kinematic evolution equation of the dynamic contact angle (1.4), we show rigorously that the solution obeys the GNBC law (2.17), if the solution has a $C^1$-regularity up to the contact line. Indeed, we show in Section 4 that the weak singularity at the contact line is removed in the GNBC model with finite $\epsilon$. We find a mesh-converging curvature at the contact line (see figure 9) and the numerical solution satisfies the GNBC law in a quasi-stationary state (i.e. for $\dot{\theta}_d = 0$). These results are consistent with the recent findings of Kulkarni et al. (2023) who demonstrated that this model indeed shows a local $C^2$-regularity at the contact line. As expected from kinematics,
the tangential stress component goes to zero at the contact line in quasi-stationary states (see figure 10). In this sense, we observe perfect apparent slip at the moving contact line. A natural follow-up of this work would be to extend the CR-GNBC to non-flat surfaces.

It is also worth noting that VOF and phase-field methods serve complementary roles in contact line modeling. Phase-field methods naturally capture diffuse interface physics at smaller scales, while VOF targets sharp-interface macroscopic simulations where computational efficiency is essential. Systematic benchmarks comparing VOF, phase-field, and molecular dynamics for sheared nanodroplets have been performed by Lācis et al. (2020, 2022), showing that both approaches can match MD results when properly calibrated, albeit through different regularization mechanisms. A direct comparison between the present CR-GNBC formulation and phase-field methods would constitute an interesting direction for future work.

We now discuss in detail the implications and scope of the specific developments achieved in this paper.

## 1. Development of the free contact angle method

We have developed a method that allows us to transport the contact angle in a kinematically consistent manner. This is a major difference to the traditional approaches of imposing constant contact angle or an angle based on a mobility law. A mobility law relates the contact angle with the contact line velocity and other fluid properties like the viscosity ratio, surface tension, surface roughness etc. Vast literature already exists on many of such mobility laws (Xia & Steen 2018; Snoeijer & Andreotti 2013b; Ludwicki et al. 2022). In steady-state wetting, where the contact angle remains constant over time, there is a natural inclination to impose a constant contact angle. At this stage, we have made the hypothesis that a constant contact angle exists at nanoscopic scales in steady-state wetting processes. The next question that arises is how to determine the value of this contact angle. Typically, this is decided by solving the Stokes flow equation while assuming a constant contact angle and predicting the interface shape as a function of the capillary number, capillary length, and the contact angle. Several well-known relations exist, such as the Cox law (Cox 1986; Voinov 1977) and the generalized Cox-Voinov law with the slip boundary condition of Chan et al. (2020). Mathematically, a simple mobility law is written as $U_{CL} = f(\theta)$, or in an inverse form $\theta = g(U_{CL})$. However, one could define generalized mobility laws such that $U_{CL} = f(\theta, \dot{\theta}, \ddot{\theta}, ...)$, or in an inverse form $\theta = g(U_{CL}, \frac{\partial u}{\partial y}, \frac{\partial^2 u}{\partial y^2}, ...)$.

Note that the generalized version of the second one involves the gradients of the velocity field at the contact line which would in turn include the outer scales. We interpret the CR-GNBC in section 2.3 as one kind of generalized mobility law. Note that unlike with a simple mobility law, we cannot impose a contact angle directly based on the contact line speed. Here, our free-angle extrapolation scheme proves beneficial. Having successfully demonstrated its capability in the CR-GNBC case, future work could explore its applicability in other forms of the generalized mobility law.

## 2. A grid independent contact region GNBC

Grid independent results are obtained for a fixed $\varepsilon$ and $\lambda$ with varying grid size $\Delta$, given that $\Delta \ll \varepsilon$ and $\Delta \ll \lambda$. We have obtained converging results even with $\varepsilon/\Delta = 5$. Obtaining grid-independent results is crucial for predicting the transition Capillary number as a function of $\varepsilon$ and $\lambda$ so that it could have a potential scope of comparison with experiments. Note that the study of Afkhami et al. (2018) was with no-slip boundary condition giving
rise to grid-dependent results in the Volume-of-Fluid framework. Such grid dependency is removed by using a Navier slip and resolving the slip length. We have shown that the curvature diverges logarithmically at the contact line for the Navier slip boundary condition. A divergence in curvature implies a divergence in the pressure field making the model physically ill-posed. Unlike the non-integrable stress singularity that results from the no-slip boundary condition, the Navier slip has an integrable singularity and hence grid-independent steady-state results can be found. A logarithmic divergence of curvature is accompanied by convergence of the contact angle in this case. With the CR-GNBC, given a fixed $\varepsilon$, we were able to confirm the finite curvature as predicted by the thin film equation 2.5 and Kulkarni et al. (2023) and get a smooth flow field. A resolution as low as five grid points in the contact region was sufficient to get converged results.

### 3. Shear stress and the GNBC law

This singularity-free behavior of CR-GNBC over the Navier slip can be seen as a result of incorporating the uncompensated Young stress. Without the uncompensated Young stress, the CR-GNBC reduces back to a classical Navier slip. From Kulkarni et al. (2023), we know that the Navier slip results in a merely continuous velocity field at the contact line. This implies that the shear stress is mathematically not defined at the contact line. Numerical results in figure 10b show that we have a non-converging spiked behavior of shear stress at the contact line. From the stream-function solution of Kulkarni et al. (2023), we can show that the shear stress remains bounded up to the contact line while the differentiability for the shear stress is lost at the contact line. We also observe, from figure 10, that once uncompensated Young stress is added, the shear stress at the contact line goes to zero in a converging and smooth manner. We numerically recover the GNBC law which relates the steady state contact angle value and the velocity at the contact line in the vanishing shear stress limit. This is perfectly in-line with requirement for a smooth flow from Fricke et al. (2018, 2019). Notably, this smooth behavior of shear stress going to zero in the CR-GNBC happens only within the $\varepsilon$ width, i.e. within the contact region. In the derivation of the GNBC from entropy principles (Fricke et al. 2020), we introduced a smoothed uncompensated Young stress as a Dirac function. This smoothed region can be seen as a physical contact region. However, for mathematical coupling of terms, it remains to be seen how the shear stress would behave if we retained a delta function GNBC formulation. That is, considering an uncompensated Young stress in the singular form of a true delta function, and how it would interact with the singularity of the Navier slip. An investigation of this case is left as a future task.

### 4. Extension to transient regime

In the paper we have dealt with the steady-state flow characteristics of the GNBC. The work should now be extended to incorporate the setups that have a transient contact line behavior. Examples of such setups include nano-scale shear droplet (Lācis et al. 2022) and a spreading drop. Karim et al. (2016) showed that despite using the same fluids and solid material, the overall configuration of the system can result in different values of the apparent contact angle even for equal contact line speeds. The forced wetting case of a plunging plate exhibited a different apparent angle than the spontaneous wetting case of a spreading drop, even when the contact line speed is equal to the plate velocity. Thus a simple mobility law cannot capture this dependence. However, the GNBC contains a term representing the rate of change of the contact angle 2.17 which makes spontaneous wetting different from the forced wetting and whether CR-GNBC could explain the experimental behavior remains an interesting question. The CR-GNBC should also be extended to setups having sustained oscillatory states. The vibrating drop is a famous example of this setup Xia & Steen (2018).
There have been several models to describe such oscillatory setups, but most of them rely on empirical relations (Kistler 1993). Sakakeeny & Ling (2021) numerically predicted the first and second modal frequencies of a vibrating droplet in two limiting cases (a) a pinned contact line and (b) a free-slip contact line. In reality, the contact line is expected to behave in between these two limits. Whether a fundamental boundary condition like the CR-GNBC could recover the modal frequencies on large scale as well as the contact line hysteresis at micrometer scale as observed by Xia & Steen (2018), remains to be tested.

# 6. Outlook

## 6.1. A non-linear generalization of the GNBC

As discussed in detail in Section 2.2, the CR-GNBC in the form

$$ -\beta(\mathbf{v}_{\parallel} - \mathbf{U}_w) = (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos\theta_d - \cos\theta_e) \mathbf{n}_{\Gamma}\delta_{\Gamma}^{\varepsilon} \quad \text{on } \partial\Omega \quad (6.1) $$

is obtained as a linear closure relation, to render the dissipation integral

$$ \mathcal{T} = \int_{\partial\Omega} ((\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos\theta_d - \cos\theta_e) \mathbf{n}_{\Gamma}\delta_{\Gamma}^{\varepsilon}) \cdot (\mathbf{v}_{\parallel} - \mathbf{U}_w) \, dA \quad (6.2) $$

non-positive. Since, according to kinematics, the viscous stress contribution vanishes in a quasi-stationary state (see Section 2), we obtain the dynamic contact angle relation

$$ -\zeta U_{cl} = \sigma(\cos\theta_d - \cos\theta_e) \quad (6.3) $$

with the contact line friction coefficient $\zeta = \beta\varepsilon$. Notably, equation (6.3) is also found in the Molecular Kinetic Theory (MKT) in the limit of low capillary number (see, e.g., Blake et al. (2015)). However, for higher capillary numbers, the MKT predicts that†

$$ U_{cl} = 2\kappa^0\Lambda \sinh[\sigma(\cos\theta_e - \cos\theta_d)/(2nk_B T)] \quad (6.4) $$

Therefore, it is interesting to formulate a closure relation for (6.2) that leads to the relation (6.4) in quasi-stationary states. Notice that (6.4) can be linearized for $U_{cl} \to 0$ using $\sinh(x) = x + O(x^3)$. Hence, the contact line friction coefficient is identified as $\zeta = (nk_B T)/(\kappa^0\Lambda)$.

For simplicity, let us assume that $U_w = 0$ in the following (the generalization to $U_w \neq 0$ is obvious). To proceed, it is useful to decompose the integral in (6.2) into its components normal and tangential to the contact line, according to

$$ \mathbf{v}_{\parallel} = (\mathbf{v}_{\parallel} \cdot \mathbf{n}_{\Gamma}) \mathbf{n}_{\Gamma} + (\mathbf{v}_{\parallel} \cdot \mathbf{t}_{\Gamma}) \mathbf{t}_{\Gamma}. $$

Here, we denote by $\mathbf{t}_{\Gamma}$ a unit tangent vector to the contact line. We obtain the representation

$$ \mathcal{T} = \mathcal{T}_{\perp} + \mathcal{T}_{\parallel} \quad (6.5) $$

with

$$ \mathcal{T}_{\parallel} = \int_{\partial\Omega} (\mathbf{t}_{\Gamma} \cdot (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel}) (\mathbf{t}_{\Gamma} \cdot \mathbf{v}_{\parallel}) \, dA $$

<footnote>† In this case, the average distance and equilibrium frequency of molecular jumps are denoted by $\Lambda$ and $\kappa^0$, respectively. Moreover, $n$ is the number of adsorption sites per unit area, $k_B$ is the Boltzmann constant and $T$ is the absolute temperature; see Blake et al. (2015) for more details.</footnote>
and

$$ \mathcal{T}_{\perp} = \int_{\partial\Omega} (\mathbf{n}_{\Gamma} \cdot (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos\theta_{d} - \cos\theta_{e}) \delta_{\Gamma}^{\epsilon}) (\mathbf{n}_{\Gamma} \cdot \mathbf{v}_{\parallel}) dA. \quad (6.6) $$

We are now looking for closure relations to ensure that $\mathcal{T}_{\perp} \le 0$ and $\mathcal{T}_{\parallel} \le 0$. A general non-linear closure for $\mathcal{T}_{\perp}$ reads as

$$ \mathbf{n}_{\Gamma} \cdot (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} + \sigma(\cos\theta_{d} - \cos\theta_{e}) \delta_{\Gamma}^{\epsilon} = -f(\mathbf{n}_{\Gamma} \cdot \mathbf{v}_{\parallel}), \quad (6.7) $$

where the scalar function $f$ satisfies the inequality

$$ xf(x) \ge 0 \quad \forall x \in \mathbb{R}. $$

Such a closure is consistent with the second law of thermodynamics because it implies that

$$ \mathcal{T}_{\perp} = - \int_{\partial\Omega} (\mathbf{n}_{\Gamma} \cdot \mathbf{v}_{\parallel}) f(\mathbf{n}_{\Gamma} \cdot \mathbf{v}_{\parallel}) \le 0. $$

The linear version of the GNBC is recovered as $f(x) = \beta x$. Analogously, a closure for the component tangential to the contact line can be obtained as

$$ \mathbf{t}_{\Gamma} \cdot (\mathbf{S}\mathbf{n}_{\partial\Omega})_{\parallel} = -f(\mathbf{t}_{\Gamma} \cdot \mathbf{v}_{\parallel}). $$

Motivated by (6.4), a special choice is

$$ f(x) = b \operatorname{arcsinh}(x/a). \quad (6.8) $$

with positive constants $a = 2\kappa^0\Lambda$ and $b = 2nk_B T$. Clearly, Equation (6.8) reduces by linearization to the original GNBC (6.1) with $\beta = b/a$ if $\mathbf{v}_{\parallel} \cdot \mathbf{n}_{\Gamma} \to 0$. Since (6.8) corresponds to the MKT (6.4) for quasi-stationary states, it may improve the standard GNBC model for higher values of the capillary number. This shall be studied in detail in the future.

## 6.2. Relaxation time scale of the GNBC

We can see from figure 7b that after $Ca_{\text{loc}}$ reaches a maximum value within a short-lived impulsive slip length time-scale $t_{\epsilon}$, there is an intermediate process of capillary relaxation before entering the slow decaying relaxation to the steady state. We refer to the timescale governing this intermediate relaxation as the GNBC relaxation time scale $t_r$. This timescale is crucial as it controls how fast a system will relax to the quasi-stationary state locally. Hence, we now examine what determines this time scale. Our formulation of the CR-GNBC model gives a contact angle evolution law 2.17 which can be re-written in terms of capillary number as

$$ \frac{\mu}{\sigma} \dot{\theta}_{d} = -\frac{Ca_{\text{cl}}}{2\lambda} + \frac{1}{2\epsilon} (\cos \theta_{d} - \cos \theta_{e}), \quad (6.9) $$

with $Ca_{\text{cl}}$ the contact line capillary number defined in the reference frame of the solid. In the lab frame of reference we can use $Ca_{\text{loc}}$ and $Ca$ to write the angle evolution as

$$ \frac{\mu}{\sigma} \dot{\theta}_{d} = \frac{Ca - Ca_{\text{loc}}}{2\lambda} + \frac{1}{2\epsilon} (\cos \theta_{d} - \cos \theta_{e}). \quad (6.10) $$

Assuming small deviations of $\theta_d$ from $\theta_e$, we get

$$ \frac{\mu}{\sigma} \dot{\theta}_{d} = \frac{Ca - Ca_{\text{loc}}}{2\lambda} + \frac{\sin \theta_{e}}{2\epsilon} (\theta_{e} - \theta_{d}). \quad (6.11) $$

Note that in above equation, the time dependent quantities as $\theta_d$, $\dot{\theta}_d$ and $Ca_{\text{loc}}$ require the solution of the Stokes flow equation subject to this boundary condition which leaves us
only with numerical tools. However, motivated by the results of our DNS, we propose an asymptotic expansion of the contact line capillary number in terms of the deviation of the contact angle away from the equilibrium angle. From this expansion, we can also see that the contact line is mainly driven by the uncompensated Young stress in this relaxation time scale. In the relaxation regime one may show that the contact-line velocity admits a regular expansion in the small angle deviation

$$ \delta = \theta_e - \theta_d, \qquad (6.12) $$

of the form

$$ \mu U_{cl} = A_0(\sigma\delta) + A_1(\sigma\delta)^2 + A_2(\sigma\delta)^3 + O(\delta^4), \qquad (6.13) $$

where $\mu$ is the dynamic viscosity, $\sigma$ the surface tension, and the coefficients $A_i$ are constants. This expansion in (6.13) applies strictly on the relaxation timescale $t_r$, over which the local capillary number $Ca_{loc}$ monotonically decreases from an initial value toward zero. Prior to this relaxation, an impulsive timescale of order $t_\epsilon$ governs the rapid rise of $Ca_{loc}$ from zero to a peak value $Ca_M$, accompanied by a corresponding change of the dynamic contact angle from $\theta_e$ to some $\theta_M$. It is important to note that (6.13) is derived in the frame in which the solid substrate is stationary. In the lab frame, where the plate translates at a constant capillary number Ca, the steady-state contact-line condition requires a finite offset in the boundary condition. The mobility law at the first order in relaxation time scale thus becomes

$$ -Ca_{loc} = A_0(\theta_e - \theta_d) + \Phi(Ca), \qquad (6.14) $$

where $\Phi(Ca)$ is the offset determined by the system's steady-state behaviour. Substituting the above expansion in (6.11), we obtain

$$ \frac{\mu}{\sigma} \dot{\theta}_d = \frac{A_0(\theta_e - \theta_d)}{2\lambda} + \frac{\sin \theta_e}{2\varepsilon}(\theta_e - \theta_d) + f(Ca), \qquad (6.15) $$

where the steady-state behaviour, or GNBC law (2.18), gives $f(Ca) = -\frac{\varepsilon}{\lambda} \frac{Ca}{\sin \theta_e}$. The solution of the above ordinary differential equation is,

$$ (\theta_e - \theta_d) = C_0 e^{-\omega t} - \frac{\varepsilon}{\lambda} \frac{Ca}{\sin \theta_e}, \qquad (6.16) $$

where

$$ \omega = \frac{\sigma}{2\mu} \left( \frac{A_0}{\lambda} + \frac{\sin \theta_e}{\varepsilon} \right). \qquad (6.17) $$

The solution tells us that the relaxation time constant $t_r = 1/\omega$ has a contribution from both the viscous stress ~$\sigma A_0/\mu\lambda$ and Young stress ~$\sigma \sin \theta_e/\mu\varepsilon$. Since $A_0$ and $\sin \theta_e$ both are of order one, the time scale is of the order of slip length-based capillary number and a contact region-based capillary number that comes out to be a few milliseconds for our flow parameters. The constant $C_0$, a function of Ca, needs to be obtained by asymptotic matching to the impulsive time response. We leave the solution at the impulsive time scale $t_\epsilon$ and a matched asymptotic solution at the full range of scales as future scope.

We now verify the existence of this relaxation timescale $t_r$ in our DNS. Re-plotting figure 7 with shifts accounting for each Ca, we see that all the curves collapse. Hence, figure 12 verifies the region's existence, which happens at the relaxation time scale, where the linear dependence in (6.14) is justified. We get $A_0 = 0.9$ from the best fit. At $t = 0$ for the relaxation time scale, we have $Ca_{loc} = Ca_M$ and $\theta_d = \theta_M$. Since the coordinates $(\theta_M, Ca_M)$ are decided by the limit of $t \to \infty$ from the impulsive time-scale, we directly use the numerical values of $\theta_M$ obtained from our simulations. Using the $C_0$ obtained from numerics for the highest
<figure><img src="image_14.png" /><figcaption>Figure 12: Collapse of shifted $Ca_{loc}$ as a function of the angle deviation $\theta_e - \theta_d$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Comparison of the dynamic angle $\theta_d$ from DNS (red curve) with the analytical solution (6.16) (black dashed curve) for the CR-GNBC case with $Ca = 0.04$, $\epsilon = 0.05$ and $\epsilon/\Delta = 40$. Inset shows the same plot in log-scale.</figcaption></figure>

resolution case of figure 8, the derived analytical solution (6.16) shows a good agreement with DNS results in figure 13.

**Author contributions.** TF, YK and MF contributed equally to this work that includes free angle method, performing simulations and writing paper with feedback from all authors. Study was performed at all three institutes hosted by SA, DB and SZ. Detailed discussions were done among all authors on all the ideas presented in the paper.

**Funding.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n° 883849). MF and DB acknowledge the financial support by the German Research Foundation (DFG) within the Collaborative Research Centre 1194 (Project-ID 265191195).

**Acknowledgments.** MF and DB gratefully acknowledge insightful discussions with Joël De Coninck (ULB).

**Declaration of interest.** The authors report no conflict of interest.

# Appendix A. Convergence study

We conduct a convergence study to demonstrate the grid independence of the CR-GNBC. In figure 14(a), we present interface shapes for a fixed $Ca = 0.12$ and $\epsilon = 0.2$ with varying resolutions, showing apparent convergence. In figure 14, we display the percentage error in the contact line position for this case, revealing second-order convergence. This confirms
<figure><img src="image_16.png" /></figure>

Figure 14: Steady-state height and interface shapes near the contact line with varying grid resolution for the CR-GNBC. In (a), $Ca = 0.12$ and $\epsilon = 0.2$ are fixed, showing convergent interface shapes. In (b), a fixed $Ca = 0.04$ reveals that due to implicit slip, steady-state solutions are achievable even with a no-slip boundary condition. Interface shapes do not converge with grid refinement, and no steady-state height is found at resolutions higher than $l_c/\Delta > 100$. (b) Percentage error in the contact line position for steady-state interface shapes obtained in (a). The reference solution is taken at 164 grid points per slip length $\epsilon/\Delta$, and the dashed lines represent second-order and first-order convergence. It is observed that above 20 grid points per $\epsilon$, a second-order convergence is achieved.

<figure><img src="image_17.png" /></figure>

Figure 15: Transition capillary number plotted against variation of (a) $\epsilon$ and $\lambda$ with $\epsilon = \lambda$, and (c) $\epsilon$ such that the slip length $\lambda$ is fixed. All simulations are carried out with the CR-GNBC and $\theta_c = 90^\circ$. The resolution for all simulations is maintained at $\min(\epsilon, \lambda)/\Delta = 5.12$.

that unlike Afkhami *et al.* (2018) our CR-GNBC method achieves grid independence for steady-state height with a fixed $Ca$ and $\epsilon$, including the $Ca_{tr}$.
**Appendix B. Transition capillary number and the contact region width ε**

Figure 15 illustrates Ca<sub>tr</sub> as a function of ε, λ and θ<sub>e</sub>. When we have ε = λ, we see that Ca<sub>tr</sub> decreases with the decrease of ε. When we break the restriction of ε = λ, we see an interesting behaviour in figure 15(b). Since we know that ε and λ, both promote smoothing behaviour, the Ca<sub>tr</sub> is decided by the larger of the two. Hence, when we decrease ε below the slip length λ, we see that Ca<sub>tr</sub> goes to a constant value. The dependence on the equilibrium contact angle shown in figure 15(b) is notably linear. While this linearity may break at smaller angles, it's important to note that our solver, which employs only horizontal heights, faces limitations in handling angles smaller than 30°.

REFERENCES

AFKHAMI, S., BUONGIORNO, J., GUION, A., POPINET, S., SAADE, Y., SCARDOVELLI, R. & ZALESKI, S. 2018 Transition in a numerical model of contact line dynamics and forced dewetting. *Journal of Computational Physics* **374**, 1061–1093.

AFKHAMI, S. & BUSSMANN, M. 2008 Height functions for applying contact angles to 2D VOF simulations. *International Journal for Numerical Methods in Fluids* **57** (4), 453–472.

AFKHAMI, S. & BUSSMANN, M. 2009 Height functions for applying contact angles to 3D VOF simulations. *International Journal for Numerical Methods in Fluids* **61** (8), 827–847.

BLAKE, T.D & HAYNES, J.M. 1969 Kinetics of liquid-liquid displacement. *Journal of Colloid and Interface Science* **30** (3), 421–423.

BLAKE, T. D. 2006 The physics of moving wetting lines. *Journal of Colloid and Interface Science* **299** (1), 1–13.

BLAKE, T. D., FERNANDEZ-TOLEDANO, J.-C., DOYEN, G. & DE CONINCK, J. 2015 Forced wetting and hydrodynamic assist. *Physics of Fluids* **27** (11), 112101.

BNÀ, S., MANSERVISI, S., SCARDOVELLI, R., YECKO, P. & ZALESKI, S. 2016 Vofi — a library to initialize the volume fraction scalar field. *Computer Physics Communications* **200**, 291–299.

BONN, D., EGGERS, J., INDEKEU, J., MEUNIER, J. & ROLLEY, E. 2009 Wetting and spreading. *Reviews of Modern Physics* **81** (2), 739–805.

CHAN, TAK SHING, KAMAL, CATHERINE, SNOEIJER, JACCO H., SPRITTLES, JAMES E. & EGGERS, JENS 2020 Cox–voinov theory with slip. *Journal of Fluid Mechanics* **900**, A8.

CHEN, XIANYANG, LU, JIACAI & TRYGGVASON, GRÉTAR 2019 Numerical simulation of self-propelled non-equal sized droplets. *Physics of Fluids* **31** (5), 052107.

COX, R. G. 1986 The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. *Journal of Fluid Mechanics* **168**, 169–194.

DEVAUCHELLE, O., JOSSERAND, C. & ZALESKI, S. 2007 Forced dewetting on porous media. *Journal of Fluid Mechanics* **574**, 343–364.

DUFFY, BRIAN R. & WILSON, STEPHEN K. 1997 A third-order differential equation arising in thin-film flows and relevant to tanner's law. *Applied Mathematics Letters* **10**, 63–68.

EGGERS, JENS 2004 Hydrodynamic theory of forced dewetting. *Phys. Rev. Lett.* **93**, 094502.

FRICKE, M. 2021 Mathematical modeling and volume-of-fluid based simulation of dynamic wetting. PhD thesis, TU Darmstadt.

FRICKE, M. & BOTHE, D. 2020 Boundary conditions for dynamic wetting - A mathematical analysis. *The European Physical Journal Special Topics* **229** (10), 1849–1865.

FRICKE, M., KÖHNE, M. & BOTHE, D. 2019 A kinematic evolution equation for the dynamic contact angle and some consequences. *Physica D: Nonlinear Phenomena* **394**, 26–43.

FRICKE, M., KÖHNE, M. & BOTHE, D. 2018 On the kinematics of contact line motion. *PAMM* **18** (1), e201800451.

FRICKE, M., MARIĆ, T. & BOTHE, D. 2020 Contact Line Advection using the geometrical Volume of Fluid Method. *Journal of Computational Physics* **407**, 109221.

FULLANA, T., ZALESKI, S. & POPINET, S. 2020 Dynamic wetting failure in curtain coating by the Volume-of-Fluid method. *The European Physical Journal Special Topics* **229** (10), 1923–1934.

FUMAGALLI, I., PAROLINI, N. & VERANI, M. 2018 On a free-surface problem with moving contact line: From variational principles to stable numerical approximations. *Journal of Computational Physics* **355**, 253–284.
DE GENNES, P.-G., BROCHARD-WYART, F. & QUÉRÉ, D. 2004 *Capillarity and Wetting Phenomena*. New York, NY: Springer New York.

GERBEAU, J.-F. & LELIÈVRE, T. 2009 Generalized Navier boundary condition and geometric conservation law for surface tension. *Computer Methods in Applied Mechanics and Engineering* **198** (5-8), 644–656.

HOCKING, L. M. 2001 Meniscus draw-up and draining. *European Journal of Applied Mathematics* **12** (3), 195–208.

HUH, C. & MASON, S. G. 1977 The steady movement of a liquid meniscus in a capillary tube. *Journal of Fluid Mechanics* **81** (03), 401–419.

HUH, C. & SCRIVEN, L. E 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. *Journal of Colloid and Interface Science* **35** (1), 85–101.

JACMIN, DAVID 2000 Contact-line dynamics of a diffuse fluid interface. *Journal of Fluid Mechanics* **402**, 57–88.

KARIM, M. A., DAVIS, S. H. & KAVEHPOUR, H. P. 2016 Forced versus spontaneous spreading of liquids. *Langmuir* **32** (40), 10153–10158, pMID: 27643428, arXiv: https://doi.org/10.1021/acs.langmuir.6b00747.

KAWAKAMI, K., KITA, Y. & YAMAMOTO, Y. 2023 Front-tracking simulation of the wetting behavior of an impinging droplet using a relaxed impermeability condition and a generalized navier boundary condition. *Computers & Fluids* **251**, 105739.

KISTLER, STEPHAN F 1993 Hydrodynamics of wetting. *Wettability* **6**, 311–430.

KULKARNI, YASH, FULLANA, TOMAS & ZALESKI, STEPHANE 2023 Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized navier boundary condition. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* **479** (2278).

LUDWICKI, JM, KERN, VR, MCCRANEY, J, BOSTWICK, JB, DANIEL, S & STEEN, PH 2022 Is contact-line mobility a material parameter? npj microgravity. *NPJ Microgravity*.

LUKYANOV, A. V. & PRYER, T. 2017 Hydrodynamics of moving contact lines: Macroscopic versus microscopic. *Langmuir* **33** (34), 8582–8590.

LĀCIS, U., JOHANSSON, P., FULLANA, T., HESS, B., AMBERG, G., BAGHERI, S. & ZALESKI, S. 2020 Steady moving contact line of water over a no-slip substrate. *The European Physical Journal Special Topics* **229** (10), 1897–1921.

LĀCIS, UĞIS, PELLEGRINO, MICHELE, SUNDIN, JOHAN, AMBERG, GUSTAV, ZALESKI, STÉPHANE, HESS, BERK & BAGHERI, SHERVIN 2022 Nanoscale sheared droplet: volume-of-fluid, phase-field and no-slip molecular dynamics. *Journal of Fluid Mechanics* **940**, A10.

MARENGO, MARCO & DE CONINCK, JOEL, ed. 2022 *The Surface Wettability Effect on Phase Change*. Springer International Publishing.

MARIĆ, T., KOTHE, D. B. & BOTHE, D. 2020 Unstructured un-split geometrical Volume-of-Fluid methods – A review. *Journal of Computational Physics* **420**, 109695.

POPINET, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. *Journal of Computational Physics* **228** (16), 5838–5866.

POPINET, S. 2015 A quadtree-adaptive multigrid solver for the serre–green–naghdi equations. *J. Comput. Phys.* **302**, 336–358.

POPINET, S. 2018 Numerical Models of Surface Tension. *Annual Review of Fluid Mechanics* **50** (1), 49–75.

POPINET, STÉPHANE & ZALESKI, STÉPHANE 1999 A front-tracking algorithm for accurate representation of surface tension. *International Journal for Numerical Methods in Fluids* **30** (6), 775–793.

PRÜSS, J. & SIMONETT, G. 2016 *Moving interfaces and quasilinear parabolic evolution equations. Monographs in Mathematics*. Switzerland: Birkhäuser.

QIAN, T., WANG, X.-P. & SHENG, P. 2003 Molecular scale contact line hydrodynamics of immiscible flows. *Physical Review E* **68** (1 Pt 2), 016306.

QIAN, T., WANG, X.-P. & SHENG, P. 2006a A variational approach to moving contact line hydrodynamics. *Journal of Fluid Mechanics* **564**, 333.

QIAN, T., WANG, X.-P. & SHENG, P. 2006b Molecular hydrodynamics of the moving contact line in two-phase immiscible flows. *Commun. Comput. Phys.* **1** (1), 1–52.

REN, W. & E, W. 2007 Boundary conditions for the moving contact line problem. *Physics of Fluids* **19** (2), 022101.

REN, W & E, W 2011a Contact line dynamics on heterogeneous surfaces. *Physics of Fluids* **23** (7), 072103.

REN, W. & E, W. 2011b Derivation of continuum models for the moving contact line problem based on thermodynamic principles. *Communications in Mathematical Sciences* **9** (2), 597–606.

SAKAKEENY, JORDAN & LING, YUE 2021 Numerical study of natural oscillations of supported drops with free
and pinned contact lines. *Physics of Fluids* **33** (6), 062109, arXiv: https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/5.0049328/15909411/062109_1_online.pdf.

SCARDOVELLI, R. & ZALESKI, S. 1999 Direct Numerical Simulation of Free-Surface and Interfacial Flow. *Annual Review of Fluid Mechanics* **31** (1), 567–603.

SHANG, XINGLONG, LUO, ZHENGYUAN, GATAPOVA, ELIZAVETA YA., KABOV, OLEG A. & BAI, BOFENG 2018 GNBC-based front-tracking method for the three-dimensional simulation of droplet motion on a solid surface. *Computers & Fluids* **172**, 181–195.

SHIKHMURZAEV, Y. D. 1993 The moving contact line on a smooth solid surface. *International Journal of Multiphase Flow* **19** (4), 589–610.

SHIKHMURZAEV, Y. D. 2006 Singularities at the moving contact line. Mathematical, physical and computational aspects. *Physica D: Nonlinear Phenomena* **217** (2), 121–133.

SHIKHMURZAEV, Y. D. 2008 *Capillary flows with forming interfaces*. Boca Raton: Chapman & Hall/CRC.

SLATTERY, J. C. 1999 *Advanced transport phenomena*. Cambridge Series in Chemical Engineering . Cambridge: Cambridge University Press.

SNOEIJER, J. H. & ANDREOTTI, B. 2013a Moving contact lines: Scales, Regimes, and Dynamical Transitions. *Annual Review of Fluid Mechanics* **45** (1), 269–292.

SNOEIJER, JACCO H. & ANDREOTTI, BRUNO 2013b Moving contact lines: Scales, regimes, and dynamical transitions. *Annual Review of Fluid Mechanics* **45** (Volume 45, 2013), 269–292.

THOMPSON, PETER A. & ROBBINS, MARK O. 1989 Simulations of contact-line motion: Slip and the dynamic contact angle. *Physical Review Letters* **63** (7), 766–769.

TRYGGVASON, G., SCARDOVELLI, R. & ZALESKI, S. 2011 *Direct numerical simulations of gas-liquid multiphase flows*. Cambridge University Press.

VOINOV, O. V. 1977 Hydrodynamics of wetting. *Fluid Dynamics* **11** (5), 714–721.

XIA, YI & STEEN, PAUL H. 2018 Moving contact-line mobility measured. *Journal of Fluid Mechanics* **841**, 767–783.

YAMAMOTO, Y., HIGASHIDA, S., TANAKA, H., WAKIMOTO, T., ITO, T. & KATOH, K. 2016 Numerical analysis of contact line dynamics passing over a single wettable defect on a wall. *Physics of Fluids* **28** (8).

YAMAMOTO, Y., ITO, T., WAKIMOTO, T. & KATOH, K. 2013 Numerical simulations of spontaneous capillary rises with very low capillary numbers using a front-tracking method combined with generalized navier boundary condition. *International Journal of Multiphase Flow* **51**, 22–32.

YAMAMOTO, Y., TOKIEDA, K., WAKIMOTO, T., ITO, T. & KATOH, K. 2014 Modeling of the dynamic wetting behavior in a capillary tube considering the macroscopic–microscopic contact angle relation and generalized navier boundary condition. *International Journal of Multiphase Flow* **59**, 106–112.

YOUNG, T. 1805 An essay on the cohesion of fluids. *Philosophical Transactions of the Royal Society of London* **95**, 65–87.

ZHANG, JIAQI & YUE, PENGTAO 2020 A level-set method for moving contact lines with contact angle hysteresis. *Journal of Computational Physics* **418**, 109636.