# Mass diffusion and bending in dynamic wetting by phase-field and sharp-interface models

Tomas Fullana, Stéphane Zaleski, Gustav Amberg

► To cite this version:

Tomas Fullana, Stéphane Zaleski, Gustav Amberg. Mass diffusion and bending in dynamic wetting by phase-field and sharp-interface models. Physical Review E, 2025, 112 (4), pp.045108. 10.1103/xmj5-dbgp. hal-05417162

HAL Id: hal-05417162

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

Submitted on 15 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mass diffusion and bending in dynamic wetting by phase-field and sharp-interface models

Tomas Fullana,¹,² Stéphane Zaleski,¹,³ and Gustav Amberg⁴,\*

¹Institut Jean le Rond d'Alembert, Sorbonne Université, CNRS, F-75005 Paris, France

²Laboratory of Fluid Mechanics and Instabilities, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland

³Institut Universitaire de France, Paris, France

⁴Flow Centre, Department of Engineering Mechanics, The Royal Institute of Technology, 100 44 Stockholm, Sweden

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

(Received 16 April 2025; accepted 21 August 2025; published 21 October 2025)

Dynamic wetting poses a well-known challenge in classical sharp-interface formulation as the no-slip wall condition leads to a contact line singularity that is typically regularized with a Navier boundary condition, often requiring empirical fitting for the slip length. On the other hand, this paradox does not appear in phase-field models as the contact line moves through diffusive mass transport. In this work, we present a toy model that accounts for mass diffusion at the contact line within a sharp-interface framework. This model is based on a theoretical relation derived from the Cahn-Hilliard equations, which links the total diffusive mass transport to the curvature at the wall. From an estimate of the chemical potential on a curved interface, we obtain an expression for the width of the highly curved region $\delta$ and the apparent angle. In the sharp-interface model, we then introduce a fictitious boundary, displaced by a distance $\delta$ into the domain, where a Navier boundary condition is applied along a dynamic apparent contact angle that accounts for the local interface bending. The robustness of the model is assessed by comparing the toy model results with standard phase-field ones on two cases: the steady state profiles of a liquid bridge between two plates moving in opposite directions and the transient behaviors of a drop spreading on a solid with a prescribed equilibrium angle. This offers a practical and efficient alternative to solve contact line problems at lower cost in a sharp-interface framework with input parameters from phase-field models.

DOI: 10.1103/xmj5-dbgp

## I. INTRODUCTION

The contact line formed at the intersection of a fluid-fluid and a solid boundary is a key element of many natural and several technological processes and its dynamics has been a subject of debate for several decades [1–7]. The classical no-slip boundary condition at a moving fluid interface leads to a force singularity at the contact line, as demonstrated by Huh and Scriven [8], implying that contact line motion would be impossible.

In the context of continuum mechanics, in sharp-interface models, one common approach is to introduce a small slip velocity at the contact line using the Navier boundary condition. The slip length, typically on the order of nanometers, is often associated with surface roughness or the mean free path of molecules [9]. The Navier slip model combined with an equilibrium contact angle on a mesh-aligned boundary has been extensively studied across various physical setups using different numerical methods. For instance, a rapidly

advancing contact line in a curtain coating scenario was simulated in a volume-of-fluid (VOF) framework [10], different numerical methods with the same contact line model were tested in a capillary rise setup [11], a withdrawing plate was simulated using the VOF method with a mesh-dependent slip length [12]. Spreading droplets were explored with mesh-dependent slip [13], and the level-set method was employed for a spreading drop case, incorporating contact line hysteresis by varying the contact angle based on the contact line speed [14]. Recently, a Generalized Navier Boundary Condition implemented within the same sharp-interface framework as the present work was shown to regularize the logarithmic singularity at the contact line by allowing a deviation of the dynamic contact angle with respect to the equilibrium one through the uncompensated Young’s stress [15,16]. These models have in common that they require empirical information on the slip length, which can be viewed as a fitting parameter.

On the other hand, in the phase-field method (PF), the Cahn-Hilliard equations are formulated from the thermodynamics of an immiscible two-component mixture and a phase function is used to represent the moving interface. The velocity field can satisfy a no-slip condition at the contact line, and the contact line moves due to the diffusive mass transport that is present in a two-component system [17,18]. While the phase-field method has the advantage of needing less empirical fitting, it is considerably more expensive. Although it has been shown that it is possible to match the results

<footnote>\*Contact author: gustava@kth.se</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by Bibsam.</footnote>
obtained on dynamic contact lines by fine-tuning the different parameters of the models [9,19], an ideal model would mimic mass diffusion by a simple model; that is to obtain results essentially equivalent to the PF ones by doing much cheaper VOF simulations.

We introduce a computational boundary at a specified short distance from the wall in the VOF simulation. There we apply a Navier condition for tangential velocity, and an apparent contact angle that accounts for the interface curvature at the wall. The basis for the model is a theoretical relation between the total diffusive mass transport across the interface and the interface curvature at the wall, that is derived from the Cahn-Hilliard equations. We compare full phase-field simulations with VOF simulations using these boundary conditions. The input parameters for these VOF simulations are the same as for the full PF solutions. The VOF solutions reproduce the PF solutions, at a much lower cost, and within the framework of well-established methods in CFD.

The paper is organized as follows. Section II describes the sharp-interface framework and Sec. III the phase-field formulation. In Sec. IV we derive the toy model for the dynamic contact line starting from the phase-field approximation of the incompressible two-phase Navier-Stokes equations. In Sec. V we show a comparison of the toy model VOF results with standard PF ones on two different systems. Conclusion follows in Sec. VI.

## II. SHARP-INTERFACE MODEL

The volume-of-fluid method for representing fluid interfaces coupled with a flow solver is well known to be suited for solving interfacial flows [20–22]. In the present work, we use the free software Basilisk [23], a platform for the solution of partial differential equations on adaptive Cartesian meshes [24–26]. We consider two incompressible, Newtonian and nonmiscible fluids. We assume that there is no mass and heat transfer at the fluid-fluid interface, that the surface tension between the two fluids $\sigma$ is constant, and neglect gravity. The two fluids will be denoted $l$ and $g$ for liquid and gas phases, with densities and dynamical viscosities $\rho_l$ and $\rho_g$, $\mu_l$ and $\mu_g$ respectively. In the one-fluid formulation, the unsteady, incompressible, variable density Navier-Stokes equations read

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \quad (1) $$

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla P + \nabla \cdot (2\mu \mathbf{D}) + \mathbf{F}_{\sigma}, \quad (2) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (3) $$

with $\mathbf{u} \equiv \mathbf{u}(\mathbf{x}, t)$ the flow velocity, $P \equiv P(\mathbf{x}, t)$ the pressure, $\mathbf{D} = (\nabla \mathbf{u} + \nabla \mathbf{u}^T)/2$ the rate of deformation tensor, $\rho \equiv \rho(\mathbf{x}, t)$ the density, $\mu \equiv \mu(\mathbf{x}, t)$ the viscosity, and $\mathbf{F}_\sigma$ the surface tension force. The surface tension term is modeled using the continuum surface force [27] so that $\mathbf{F}_\sigma = \sigma \kappa \nabla \chi$ with $\kappa$ the curvature, $\mathbf{n}$ the normal unit vector to the interface and $\chi$ the Heaviside function

$$ \chi(\mathbf{x}, t) = \begin{cases} 0, & \text{in the gas phase,} \\ 1, & \text{in the liquid phase,} \end{cases} \quad (4) $$

also used to locate the interface and define the viscosity and density [22]. The volume fraction $c(\mathbf{x}, t)$ corresponding to the average value of $\chi$ in each cell is used to replace the advection of the density (1) by the following advection equation:

$$ \frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0. \quad (5) $$

The local density and viscosity are then defined from the local volume fraction $c$

$$ \rho(c) \equiv c \rho_l + (1-c) \rho_g, \quad (6) $$

$$ \mu(c) \equiv c \mu_l + (1-c) \mu_g. $$

The incompressible Navier-Stokes equations are approximated using a time-staggered approximate projection method on a Cartesian grid. The advection term is discretized with the explicit and conservative Bell–Colella–Glaz second-order unsplit upwind scheme [28]. For the discretization of the viscous diffusion term, a second-order Crank–Nicholson fully-implicit scheme is used. Spatial discretization is achieved using a quadtree (in two dimensions) adaptive mesh refinement on collocated grids. The adaptive mesh technique is based on the wavelet decomposition of the variables such as the velocity field or the volume fraction and is efficient and critical to the success of the simulations as it allows high resolution only where needed, therefore reducing the cost of calculation.

All the variables are colocated at the center of each square discretization volume. Consistently with a finite-volume formulation, the variables are interpreted as the volume-averaged values for the corresponding discretization volume. To solve the advection equation the geometrical VOF scheme is used and proceeds in two steps: (i) Interface reconstruction and (ii) Geometrical of flux estimation and interface advection. The reconstruction is a ‘piecewise linear interface calculation’ (PLIC), followed by a Lagrangian advection.

The resolution of the surface tension term is directly dependent on the accuracy of the curvature calculation. The Height-Function methodology is a VOF-based technique for calculating interface normals and curvatures [12,29]. 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 3×7 stencil around an interface cell is constructed and the heights are evaluated by summing volume fractions vertically

$$ h_j = \sum_{k=i-3}^{k=i+3} c_{j,k} \Delta v, \quad (7) $$

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

$$ \mathbf{n} = (h_x, -1), \quad (8) \\ \kappa = \frac{h_{xx}}{(1+h_x^2)^{3/2}}, $$

where $h_x$ and $h_{xx}$ are discretized using second-order central differences. This method is also used to model the contact angle at the triple point. 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 VOF 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. The boundary conditions in the VOF toy model, derived in Sec. IV, will result in a simple slip model coupled to a dynamic contact angle. The slip model, or Navier boundary condition, has been implemented and used in various studies within this numerical framework [9,10]. In two dimensions, the Navier boundary condition reads

$$ 
\begin{aligned}
u_x|_{y=0} + \lambda \left. \frac{\partial u_x}{\partial y} \right|_{y=0} &= U, \\
u_y|_{y=0} &= 0,
\end{aligned} 
\quad (9) $$

with $u_x|_{y=0}$ and $u_y|_{y=0}$ the x and y component of the velocity at the solid boundary, $\lambda$ the slip length, and $U$ the prescribed velocity of the moving solid. Note that (9) is written with the boundary located at $y = 0$; the location of the computational boundary will be lifted in the toy model (Sec. IV) such that the boundary is located at $y = +\lambda$. Moreover, it is important to remark that the accuracy of the numerical method is dependent on the resolution of the smallest length scale, the slip length $\lambda$. In the following numerical setups, the grid size will be chosen such that $\lambda/\Delta v \simeq 4$. This particular choice is based on our experience in contact line systems with the VOF method where sufficient resolution (i.e., grid size smaller than the slip length) is necessary to achieve grid-converged results [10].

### III. DIFFUSE INTERFACE MODEL

In the phase-field method, the governing equations are derived from the thermodynamic potentials of the system, together with the assumption of a surface energy associated with an interface [18,30,31]. It is, therefore, possible to consider different physical situations with relative ease. It is also straightforward to implement numerically, since interfaces are not tracked explicitly. Instead, a variable is introduced that has different constant values in the two phases with a steep transition between the two in the diffuse interface [31,32]. One of the major drawbacks of this method, however, is that the width of the interface $\epsilon$ must be small to match the proper interface dynamics.

Similarly to the one-fluid formulation, the phase-field model introduces a phase variable $C(\mathbf{x}, t)$ ranging from 1 to -1. This phase-field variable $C(\mathbf{x}, t)$ is governed by a convection-diffusion equation

$$ \frac{\partial C}{\partial t} = \nabla \cdot (-\mathbf{F}_d - \mathbf{F}_c), \quad (10) $$

where $\mathbf{F}_d$ is the diffusive flux and $\mathbf{F}_c$ is the convective flux. The latter, in the context of incompressible flows, takes the simple form

$$ \mathbf{F}_c = \mathbf{u} \cdot \nabla C, \quad (11) $$

and the diffusive flux is

$$ \mathbf{F}_d = -M \nabla \phi, \quad (12) $$

where $M$ is a proportionality coefficient called the phase-field mobility and $\phi$ the chemical potential defined as

$$ \phi = \beta \Psi'(C) - \alpha \nabla^2 C. \quad (13) $$

In the chemical potential, we have two parameters $\alpha$ and $\beta$, which are related to the surface tension $\sigma$ and the characteristic thickness of the diffuse interface $\epsilon$ as

$$ \sigma = \frac{2}{3} \sqrt{2\alpha\beta}, \quad \epsilon = \sqrt{\frac{\alpha}{\beta}}. \quad (14) $$

In addition, it contains the derivative of the standard double-well potential

$$ \Psi(C) = \frac{(C+1)^2(C-1)^2}{4}. \quad (15) $$

With the definitions above (10)–(15), the convection-diffusion equation, referred to as the Cahn-Hilliard [30] equation, can be written as

$$ \frac{\partial C}{\partial t} = \nabla \cdot [M \nabla(\beta \Psi'(C) - \alpha \nabla^2 C)] - \mathbf{u} \cdot \nabla C. \quad (16) $$

Equation (16) is the phase-field analog of (5). The boundary conditions for this fourth-order partial differential equation will be detailed in our study of dynamic wetting by the phase-field method. Taking the same notation as in (2), the momentum equation in the Navier-Stokes equations becomes

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla P + \nabla \cdot (2\mu \mathbf{D}) + \phi \nabla C, \quad (17) $$

where $\phi \nabla C$ corresponds to the surface tension force and acts over the diffuse interface region. Moreover, the density and viscosity are now defined through the phase-field variable $C$ as

$$ \begin{aligned} \rho(C) &\equiv \rho_l \frac{C+1}{2} + \rho_g \frac{C-1}{2}, \\ \mu(C) &\equiv \mu_l \frac{C+1}{2} + \mu_g \frac{C-1}{2}. \end{aligned} \quad (18) $$

We now state the boundary conditions in the presence of a contact line on a solid, characterized by its equilibrium angle (the angle between the fluid-fluid interface and the wall). The convection-diffusion equation (16) is a fourth-order partial differential equation and requires two boundary conditions. First, we impose a nonequilibrium wetting condition [31,33] on the solid wall,

$$ -\mu_f \epsilon \left( \frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C \right) = \alpha \nabla C \cdot \mathbf{n} - \sigma \cos \theta_e g'(C), \quad (19) $$

where $\mu_f$ is a contact line friction parameter, having the same units as bulk dynamic viscosity. Here, $\theta_e$ is the equilibrium contact angle and

$$ g(C) = \frac{1}{4}(2 - 3C + C^3), \quad (20) $$

is a switching function describing a smooth transition between both phases. The unit normal vector $\mathbf{n}$ is directed from the fluid to the surrounding solid.

If one sets $\mu_f = 0$, the contact angle is always enforced to the equilibrium angle $\theta_e$. This will be the case in the toy model when comparing PF to VOF simulations. Nonzero contact line friction allows the dynamic contact angle to evolve naturally
as a function of contact line speed. The second boundary condition for the phase function is a zero diffusive flux of chemical potential through the boundaries

$$ \nabla \phi \cdot \mathbf{n} = 0. \qquad (21) $$

The phase-field equations are written in a weak form, and solved using the open-source finite-element software FreeFEM by Hecht, Pironneau, and others [34,35], allowing the specification of a finite-element weak form as input. Linear elements were used for the PF variables, while the fluid flow was resolved using Taylor-Hood elements (quadratic for velocity and linear for pressure). The numerical code is available in the Github repository FreeFEM-NS-CH [36]. Similarly to the VOF implementation, the accuracy of the method depends on the resolution of the smallest length scale of the system, the interface width $\epsilon$. In the PF setups, the grid size will be chosen such that $\epsilon/\Delta_p \simeq 16$.

## IV. THE TOY MODEL: RELATING MASS FLUX TO CURVATURE

We start by formulating the dimensionless incompressible two-phase Navier-Stokes equations in the phase-field approximation. As reference quantities we choose $U, L, L/U$, and $\mu_l U/L$ for velocity, length, time, and pressure, with $L$ denoting a macroscopic length, for example the droplet radius, and $U$ the contact line speed. The dimensionless numbers that appear are the Reynolds number $Re = \rho_l U L / \mu_l$, the Peclet number $Pe = U L / D$, where $D = M \sigma / \epsilon$ is the mass diffusion coefficient, and the capillary number $Ca = \mu_l U / \sigma$. Lastly, the Cahn number—or dimensionless interface width—is defined as $Cn = \epsilon / L$ and $Ca_F = \mu_F U / \sigma$ is the line friction capillary number, based on the line friction parameter $\mu_F$. In the rest of this study, we consider the line friction coefficient to be zero ($\mu_F = 0$), in order to enforce the dynamic contact angle to the equilibrium one. The dimensionless equations therefore read

$$ \rho_l \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\frac{1}{\mathrm{Re}} \nabla P^* + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u} - \frac{1}{\mathrm{Re} \mathrm{Cn} \mathrm{Ca}} C \nabla \phi, \qquad (22) $$

$$ \frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C = \frac{1}{\mathrm{Pe}} \nabla^2 \phi, \qquad (23) $$

with

$$ \phi = -\mathrm{Cn}^2 \nabla^2 C + \Psi'(C), \qquad (24) $$

and

$$ \Psi'(C) = (C^2 - 1)C. \qquad (25) $$

Here the pressure has been replaced by the reduced pressure $P^* = P - C\phi/(\mathrm{Cn}\mathrm{Ca})$, which is continuous across the interface at equilibrium.

The wetting boundary condition, with $\mu_F = 0$, is now

$$ \frac{3}{2\sqrt{2}} \mathrm{Cn} \nabla C \cdot \mathbf{n} = \cos(\theta_e) g'(C). \qquad (26) $$

In what follows, we assume that there is a small region near the contact line where the interface has locally a high curvature, as shown in Fig. 1. We assume the contact angle at

<figure><img src="image_3.png" /><figcaption>FIG. 1. Sketch of the region of high curvature close to a moving contact line.</figcaption></figure>

the wall to be kept at the nominal equilibrium value $\theta_e$, but that it will adjust to an apparent contact angle $\theta_a$ over a short distance $\delta$. We will now proceed to derive estimates for the chemical potential on a curved interface, and from this, the local curvature at the moving contact line. Using this, we will arrive at semiquantitative expressions for the width of the highly curved region, and the value of the apparent contact angle $\theta_a$. Assume a length scale $\ell$, a viscous balance in the momentum equation, and that the velocity is properly scaled (dimensionless $|\mathbf{u}| \sim 1$). Estimates of the relevant terms in (22)–(23) then give

$$ \frac{1}{\mathrm{Re}} \frac{1}{\ell^2} \sim \frac{1}{\mathrm{Ca} \mathrm{Cn} \mathrm{Re}} \frac{\phi}{\ell}, \qquad (27) $$

and

$$ \frac{1}{\ell} \sim \frac{1}{\mathrm{Pe}} \frac{\phi}{\ell^2}. \qquad (28) $$

Solving for $\ell$ yields

$$ \ell \sim \sqrt{\frac{\mathrm{Ca} \mathrm{Cn}}{\mathrm{Pe}}} = \frac{\sqrt{\mu_l M}}{L}, \qquad (29) $$

which agrees with the scaling mentioned by several authors previously [18,31,37]. Briant et al. [38] have given a somewhat different expression $(\mu_l M \epsilon^2)^{1/4}/L$. This however involves the interface thickness in the estimates of the concentration gradients, while the length scale in (29) characterizes the length scale for diffusion in the contact line region, away from the interface.

We now derive an expression for the chemical potential on a locally curved interface. Assuming a planar geometry and a circular interface of radius $R$, the chemical potential as defined in (24) can be written as

$$ \phi = -\mathrm{Cn}^2 \left( \frac{\partial^2 C}{\partial r^2} + \frac{1}{r} \frac{\partial C}{\partial r} \right) + \Psi'(C). \qquad (30) $$

$C$ is expected to approach +1 inside the droplet, and -1 outside. Since the dimensionless interface width is $\mathrm{Cn} \ll 1$, the solution to this equation has a boundary layer character, and we introduce a new radial coordinate $s$ such that $r = R + \mathrm{Cn} \cdot s$ and expand (30) at zeroth and first order

$$ C = C_0 + \mathrm{Cn} \cdot C_1 + O(\mathrm{Cn}^2), \qquad (31) $$

$$ \phi = \phi_0 + \mathrm{Cn} \cdot \phi_1 + O(\mathrm{Cn}^2). \qquad (32) $$
The problems for the zeroth and first order are now

$$ \phi_0 = -\frac{\partial^2 C_0}{\partial s^2} + \Psi'(C_0), \quad (33) $$

and

$$ \phi_1 = - \left( \frac{\partial^2 C_1}{\partial s^2} - \frac{1}{R} \frac{\partial C_0}{\partial s} \right) + C_1 \Psi''(C_0). \quad (34) $$

The zeroth order (33) gives the equilibrium solution as

$$ C_0 = -\tanh(s/\sqrt{2}), \quad \text{with } \phi_0 = 0. \quad (35) $$

Multiplying (34) by $\partial C_0/\partial s$ and partially integrating, a solvability condition is obtained for the first-order value of the chemical potential

$$ \int \phi_1 C_{0,s} ds = - \int (C_{0,s} C_{1,ss} + C_1 (\Psi'(C_0))_{,s}) ds \\ - \frac{1}{R} \int (C_{0,s})^2 ds, \quad (36) $$

where the subscript denotes derivative with respect to s. The chemical potential is continuous at the interface, so the first order chemical potential $\phi_1$ is a constant in (34) and (36). By partially integrating the first term on the right hand side that term is seen to be zero. Using that

$$ \int C_{0,s} ds = -2 \quad \text{and} \quad \int (C_{0,s})^2 ds = \frac{2\sqrt{2}}{3}, \quad (37) $$

we obtain the value of the chemical potential at the interface $\phi_i$

$$ \phi_i = Cn \phi_1 = \frac{1}{2} \frac{2\sqrt{2}}{3} \frac{Cn}{R}. \quad (38) $$

Turning now to (23), which expresses the mass balance of the concentration C, the solution to this equation varies on two different length scales: one is the interface width and the other a larger length scale characterizing the region around the contact line where mass diffusion restores the concentration to its equilibrium value, i.e., the Yue-Zhou-Feng length-scale $\ell$ in (29). In order to establish a quantitative model, we now rewrite the scaling of (23) in (28) as an equality

$$ b \frac{2}{\delta} = \frac{1}{Pe} \frac{\phi_i}{\delta^2}. \quad (39) $$

Here we introduced a nondimensional factor b of order unity, which accounts for the shape of the velocity profile at the contact line. We note that the full Eq. (23), when rescaled according to (29), will be a nondimensional equation, independent of model parameters. As b is determined from a solution of this rescaled system, we conclude that b is independent of the model parameters; with the exception of the equilibrium contact angle $\theta_e$. For partially wetting systems where $\theta_e$ is not close to 0 or 180° we expect this dependence to be weak.

We now consider the Navier-Stokes equations (22). Assuming a balance between the viscous term and the capillary term, we estimate the relevant terms in the Navier-Stokes equation as

$$ \frac{1}{\text{Re}} \frac{\tau}{\delta} = \frac{1}{\text{Ca Cn Re}} \phi_i \frac{2}{\delta}. \quad (40) $$

Here a dimensionless parameter a of order unity is introduced to account for the viscous stress $\tau = au/\delta$. As discussed above in relation to the nondimensional parameter b, a is similarly determined from a rescaled equation, free from model parameters, as long as the flow is considered inertialess. However, if the two fluid regions have different viscosities, the viscosity ratio does enter. From our numerical experiments however this dependence appears weak.

Solving now for R and $\delta$ from (38)–(40) yields the local curvature as a function of the flow parameters

$$ R = \frac{1}{2} \frac{2\sqrt{2}}{3} \frac{1}{\sqrt{ab}} \sqrt{\frac{Cn}{PeCa}}. \quad (41) $$

$$ \delta = \frac{1}{2} \sqrt{\frac{a}{b}} \sqrt{\frac{CaCn}{Pe}}. \quad (42) $$

Throughout the treatment above, the viscosity is taken as the higher liquid viscosity, in the case of unequal viscosities in the two regions, and as noted above, the dependence of our results on the viscosity ratio is weak. This may seem to be in contrast with the observation of Yue et al. [18] (Figs. 9 and 11), who note that the distance from the wall to the stagnation point near the contact line (the boundary layer thickness) should be based on the geometric mean viscosity $\mu_e = \sqrt{\mu_l \mu_a}$. Here, however, the interpretation of the boundary layer thickness is not the distance to the stagnation point per se, but rather estimates the width of the highly curved region. One explanation could be that, as is clear from flow fields in the phase-field simulations, the stagnation point is displaced into the low viscosity region when the viscosities differ. For the analysis in the present work, it is the action of the viscous stresses and the shape of the velocity profile in the contact line region that is important, and not the location of the stagnation point as such.

The local radius of curvature R at the contact line allows us to estimate how the local angle $\theta_a$ varies as a function of the distance to the wall. Assuming a locally circular interface of radius R according to (41), some geometrical considerations show that the angle $\theta_a$, taken at a distance $\delta$ from the wall according to (42), is given by the following relation:

$$ \cos(\theta_a) = \cos(\theta_e) + \frac{\delta}{R} = \cos(\theta_e) + \frac{3}{2\sqrt{2}} \text{Ca} a. \quad (43) $$

In what follows we will use (43) to calculate the apparent contact angle $\theta_a$ from the local capillary number Ca, and apply this at a fictitious boundary displaced by the distance $\delta$ into the flow domain, as indicated in Fig. 1.

# V. RESULTS

## A. Sheared drop

In the first test case, we study the steady shape of a nanoscopic droplet confined between two moving walls similarly to Lacis et al. [9,19]. Initially, we use the phase-field model to solve the system in order to validate the correlation between the radius of curvature and the flow parameters, as expressed in (41). We choose a constant capillary number $Ca = 0.0212$, which is below the critical capillary number, while systematically varying the Peclet and Cahn numbers. The equilibrium angle is set to $\theta_e = 90^\circ$ for both the
<figure><img src="image_4.png" /><figcaption>FIG. 2. Correlation between the radius of curvature R and (41). (a) Each simulation corresponds to a unique set of Peclet (Pe) and Cahn (Cn) numbers while maintaining a constant capillary number (Ca). The toy model parameters are set to a = 2.5 and b = 0.2170, and the fitted line demonstrates a clear one-to-one relationship between these two quantities. (b) Steady-state solutions of phase-field (PF) simulations for cases highlighted in red in the left panel.</figcaption></figure>

advancing and receding contact lines, and the geometrical factor a = 2.5 is kept fixed. For each simulation run, we compute the radius of curvature R once the system reaches steady state, using the mean curvature at the contact line. The results, illustrated in Fig. 2, demonstrate a remarkable agreement between R and (41), establishing a clear one-to-one relationship between these two quantities. In this preliminary analysis, the parameter b = 0.2170 is fixed (see Appendix B for details). In Sec. V B, we will further substantiate the appropriateness of this choice in the context of a spreading droplet.

We now test our toy model on this setup. Using the relations described in Sec. IV, we compare full PF simulations—in the sense that the full domain is solved—with VOF simulations using the toy model approximation. In the toy model, the bottom boundary is now a computational boundary that is located a distance +δ away from the wall. The domain comprised between the wall y = 0 and the computational boundary y = δ is not simulated. There, we apply the dynamic angle relation (43), where the apparent angle θ_a is now a dynamic contact angle θ_d

$$ \theta_d = \cos^{-1} \left( \cos \theta_e + \frac{3}{2\sqrt{2}} Ca a \right), \quad (44) $$

with Ca the contact line capillary number. Moreover, in order to simulate a true no-slip condition at the wall, we use δ as the slip length for the tangential velocity (again located at the

computational boundary)

$$ u_x|_{y=\delta} - \delta \frac{\partial u_x}{\partial y} \Big|_{y=\delta} = \pm U_w, \quad (45) $$

with U_w the wall velocity. Figure 3 shows the resulting steady state interfaces for the left side of a sheared droplet, for a range of Pe and Cn. The full PF simulations are compared to the corresponding VOF results, using Eqs. (44) and (45) as boundary conditions, with good agreement. A more discerning test is shown in Fig. 4 where the local angle that the interface makes with the horizontal is shown as a function of vertical position, for the same six cases as in Fig. 3. The details of the profile shapes are near identical in all cases.

## B. Spreading drop

In this second test case, we study a simple drop spreading over a no-slip substrate. We impose an equilibrium angle which governs the final shape of the drop at equilibrium. Contrary to the previous case, we need to relax the impermeability condition by allowing a controlled mass flux through the computational boundary. This introduced vertical sinking velocity, denoted as v_s, is linked to the net mass change, being either a loss (in the case of θ_e < θ_0) or a gain (in the case of θ_e > θ_0) of mass, with θ_0 the initial angle of the drop. The boundary condition governing the normal component of
<figure><img src="image_5.png" /><figcaption>FIG. 3. Comparison between the phase-field left interface (black points) and the volume-of-fluid interface (red line) for six distinct sets of Peclet (Pe) and Cahn (Cn) numbers. Panels (a)–(f) correspond to the following parameter sets: (Pe = 30, Cn = 0.02), (Pe = 30, Cn = 0.04), (Pe = 100, Cn = 0.02), (Pe = 160, Cn = 0.04), (Pe = 300, Cn = 0.02), and (Pe = 300, Cn = 0.04).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4. Comparison between the phase-field angle along the left interface (black points) and the volume-of-fluid angle (red symbols) for six distinct sets of Peclet (Pe) and Cahn (Cn) numbers.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Comparison of interface shapes at $t = 0.1$ for phase-field simulations (black lines) and volume-of-fluid toy model (red lines). (a)–(c): Effect of the variation of the geometrical factor $a$ for a fixed value of $b$. (d)–(e) Effect of the variation of $b$ for a fixed value of $a$.</figcaption></figure>

velocity at the wall is now

$$u_y|_{y=\delta} = v_s. \qquad (46)$$

Details on the numerical implementation of $v_s$ are provided in Appendix A.

The initial radius of the drop is $R_0 = 0.5$ with an initial contact angle $\theta_0 = 90^\circ$. The equilibrium angle is set to $\theta_e = 70^\circ$. The dimensionless parameters governing the flow are $Ca = 0.0212$, $Re = 3.978$, $Cn = 0.01$, and $Pe = 1$. Figure 5 shows the influence of varying the parameters $a$ and $b$ for this transient case. Each panel shows a close-up of the left contact line region, with PF results for three different times, and the corresponding VOF results, using a range of values for $a$ and $b$. The best match is found for $a = 3$ and $b = 0.2170$ confirming the values chosen in the previous case. Note that the difference in $a$ between both cases is attributed to the different equilibrium contact angle. Figure 6 shows the diameter of the drop as a function of time for the PF and the VOF simulations. In order to compare the displacement, the position of the contact line in the VOF setup is projected back to the wall $y = 0$ using the extracted apparent angle located at $y = \delta + \Delta$, in order to avoid the numerical difference between imposed and extracted angles $\Delta\theta_{num}$ at $y = \delta$. From the comparison between both models, we can observe that the displacement curves superimpose almost perfectly. The final equilibrium position of the contact line is the same in both cases, implying that the sinking velocity $v_{sink}$ allowed

us to accurately represent the mass flux through the computational boundary. In the transient state, however, we see a slight difference between times 0.1 and 0.3 where the VOF simulation overshoots the solution. This behavior is similar to ones observed previously in slip models [39].

Figure 7 shows both the PF and VOF interfaces at $t = 0$, 0.05, 0.1, 0.15, 0.2, 0.4, and Fig. 8 shows a close-up of the same interfaces near the contact line. The solid blue line

<figure><img src="image_8.png" /><figcaption>FIG. 6. Dimensionless diameter as a function of time for the phase-field and toy model volume-of-fluid simulations for the best values of $a$ and $b$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Phase-field and toy model volume-of-fluid interfaces at times $t = 0, 0.05, 0.1, 0.15, 0.2, 0.4$. The blue line corresponds to the location of the computational boundary in the toy model.</figcaption></figure>

indicates the VOF computational boundary: as in the previous case, the domain is truncated at $y = \delta$, so the region $0 < y < \delta$ is not simulated. The superposition of the PF and VOF interfaces demonstrates overall good agreement in position and shape. A closer look in Fig. 8 reveals a discontinuity in the final VOF segment close to the contact line. This feature manifests the numerical discrepancy $\Delta\theta_{\text{num}}$, the difference between the contact angle actually reconstructed by the VOF procedure and the angle imposed via the height functions in the ghost layer.

## VI. CONCLUSION

The phase-field model for dynamic wetting simulations has the advantage that it can be explicitly connected to a particular physical situation, i.e., that of a mixture of two near immiscible fluids, and thus in principle only depends on physically meaningful parameters. Computationally it is however quite expensive, and it is not suited for large scale

<figure><img src="image_10.png" /><figcaption>FIG. 8. Enlarged view around the contact line of the phase-field and toy model volume-of-fluid interfaces at times $t = 0, 0.05, 0.1, 0.15, 0.2, 0.4$. The blue line corresponds to the location of the computational boundary in the toy model.</figcaption></figure>

CFD applications. On the other hand, the volume-of-fluid method which is the method of choice for such simulations, depends on *ad hoc* or numerical parameters for contact line motion.

In this paper we have attempted to bridge this gap by making a local analysis of the phase-field equations near the contact line, establishing the connection between a local strong curvature and the contact line speed. VOF simulations are generally successful at capturing the interface in the bulk of the liquid, and it is only in the region near the contact line that the special features of the PF are important. Using this local analysis we have parametrized the PF results, and put them in the form of a boundary condition that can be used in VOF simulations of moving contact lines. We thus propose that this will allow VOF simulations to give results that are equivalent in many respects to the results of a full PF simulation, but at the cost of a VOF simulation. We note in Appendix C that the common case of a single component phase-change fluid can be treated in a similar way, and that an equivalent model for the contact angle can be derived, see Eq. (C17).

## ACKNOWLEDGMENTS

S.Z. and T.F. thank the FLOW laboratory of KTH for hosting them and Shervin Bagheri for funding the visits in 2019, which allowed this collaboration to start. Stephane Zaleski was funded by the European Research Council (ERC) through the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 883849 TRUFLOW).

## DATA AVAILABILITY

The data that support the findings of this article are openly available [40].

## APPENDIX A: SINK VELOCITY

In the toy model, a computational boundary is introduced at a distance $+\delta$ from the wall where we apply a Navier boundary condition. When considering a spreading drop, the impermeability condition ($u_y|_{y=\delta} = 0$) needs to be relaxed in order to take into account the mass flux through that computational boundary. To that end, we introduce a sink velocity $v_{\text{sink}}$, such that

$$u_y|_{y=\delta} = v_{\text{sink}}.$$

To calculate the amount of mass lost during one time step, we compute the area of the polygon formed by the two contact line points, at times $t^{n-1}$ and $t^n$, and their linear projection onto the actual wall (at $y = 0$). In practice, the angle used for the projection is an apparent angle located one cell above the contact line. The computed area is then translated to the sink velocity by the simple formula

$$v_{\text{sink}} = \frac{1}{\tau} \frac{A}{L},$$

where $\tau$ is the time step, $L$ the length of the wall, and $A$ the signed area—defined later. Algorithm 1 summarizes the required steps at a given instant.
ALGORITHM 1. Sink velocity pseudocode.

for each cell do

Locate the cell one grid point above the contact line

Compute the contact line position $x_{CL}^n$ at $y = \delta + \Delta$

Compute the apparent angle $\theta_{app}$ using the unit normal $\mathbf{n}$

Project the contact line position

$$\tilde{x}_{CL}^n = x_{CL}^n + (\delta + \Delta) \tan(\frac{\pi}{2} - \theta_{app})$$

Compute the polygon area $A$ formed by the points

$$(x_{CL}^n, \delta + \Delta), (\tilde{x}_{CL}^n, 0), (\tilde{x}_{CL}^{n-1}, 0), (x_{CL}^{n-1}, \delta + \Delta)$$

where $x_{CL}^{n-1}$ and $\tilde{x}_{CL}^{n-1}$ were stored in the previous time step
Determine the sink velocity

$$v_{sink} = \frac{1}{\tau} \frac{A}{L}$$

Store $x_{CL}^{n-1} = x_{CL}^n$ and $\tilde{x}_{CL}^{n-1} = \tilde{x}_{CL}^n$

end

We validate the method by considering a drop spreading over a substrate for different values of $\delta$. The drop is initialized with a diameter $D = 1$ in a 2×2 domain and we apply a constant contact angle $\theta_e = 70^{\circ}$. The viscosities, densities, and the surface tension are set to 1. The simulations are run until the equilibrium position is reached. In Fig. 9, we show the total mass as a function of time when applying the sink velocity for $\delta = 0.2, 0.05, 0.02$. The initial mass inside the drop decreases as spreading occurs. In Fig. 10, we compare the final interfaces for three cases: (i) toy model with sink velocity, (ii) toy model without sink velocity, and (iii) full domain. From the results, we can observe the sink velocity allows us to match to the final true equilibrium shape (full domain). The error in final interface is a function of $\delta$. As $\delta$ is increased, the linear extrapolation procedure—that does not take into account the curved interface at the contact line—induces a higher difference in final shapes. In practice, when using the toy model on a real setup, $\delta$ is small compared to the domain size.

<figure><img src="image_11.png" /><figcaption>FIG. 9. Total mass of the drop as a function of time with or without sink velocity for different computational boundaries placed at $+\delta$ with $\delta = 0.02, 0.05, 0.2$.</figcaption></figure>

APPENDIX B: TOY MODEL CALCULATION OF PARAMETER b

Rescaling the steady version of (23) using a length-scale $\ell$ and a magnitude $\phi \sim \text{Pe} \ell$ according to (28), we get

$$\mathbf{u} \cdot \tilde{\nabla} C = \tilde{\nabla}^2 \tilde{\phi}, \qquad (B1)$$

where $\tilde{\nabla} = \ell \nabla$ and $\tilde{\phi} = \phi/(\text{Pe} \ell)$. This equation is thus universal in the sense that all dependencies of Ca, Pe, Cn, etc. are scaled out and Eq. (B1) could be solved once and for all. The only remaining independent input parameter is the equilibrium contact angle $\theta_e$. The correspondingly rescaled version of the momentum equation will introduce also the ratio of the viscosities in the two fluid regions.

In order to arrive at a simple quantitative model, the dimensionless order unity parameter $b$ is introduced in Eq. (39). This parameter is expected to essentially capture the equivalent of solving Eq. (B1), and we thus expect that it may depend on $\theta_e$ and the viscosity ratio, but that it is independent of Ca, Pe, Cn, etc.

As a bold attempt at estimating $b$ the following model equation was analyzed:

$$\text{Pe} \delta \frac{dC}{dx} = \frac{d^2}{dx^2} \left( -\frac{d^2C}{dx^2} - \frac{Cn}{R} \frac{dC}{dx} + \Psi'(C) \right). \qquad (B2)$$

Similarly to the treatment in (30)–(37), we consider a perturbation solution around the equilibrium solution in (35). The equation for the first-order perturbation becomes

$$-\frac{\text{Pe} \delta}{\text{Cn}} \frac{dC_0}{dx} + \frac{1}{R} \frac{d^3 C_0}{dx^3} = \frac{d^2}{dx^2} \left( -\frac{d^2 C_1}{dx^2} + \Psi''(C_0) C_1 \right). \qquad (B3)$$

Multiplying both sides by an $F(x)$ such that $\frac{dF(x)}{dx} = C_0(x)$, and repeatedly integrating the right hand side partially, it is seen that the right hand side vanishes, yielding a solvability condition on the remaining terms:

$$\frac{\text{Pe} \delta}{\text{Cn}} \int F \frac{d^2 F}{dx^2} dx = \frac{1}{R} \int F \frac{d^4 F}{dx^4} dx. \qquad (B4)$$

By rearranging this, using (39), we obtain

$$\frac{\text{Pe} \delta}{\text{Cn}/R} = \frac{1}{2b} \frac{1}{2} \frac{2\sqrt{2}}{3} = \frac{\int F \frac{d^4 F}{dx^4} dx}{\int F \frac{d^2 F}{dx^2} dx}. \qquad (B5)$$

The integrals on the right hand side can both be calculated, giving $\int F \frac{d^4 F}{dx^4} dx = \frac{2\sqrt{2}}{3}$ and $\int F \frac{d^2 F}{dx^2} dx = 2\sqrt{2}(1 - \ln 2)$.

This yields

$$b = \frac{1 - \ln 2}{\sqrt{2}} = 0.2170. \qquad (B6)$$

It should be emphasized that the derivation of this number is by no means a rigorous proof, and that $b$, as well as the parameter $a$ introduced in (40), are to be adjusted according to comparisons between the model in (41)–(43), and full solutions of the phase field equations, as shown in Fig. 5. As seen there however, the value $b = 0.2170$ gives satisfactory results, and it has been used throughout.
<figure><img src="image_12.png" /><figcaption>FIG. 10. Validation of the sink velocity imposed at the computational boundary placed at +δ with δ = 0.02, 0.05, 0.2 denoted by the dotted black lines in the figures. The equilibrium angle is θ$_{e}$ = 70° for the three cases considered.</figcaption></figure>

## APPENDIX C: SINGLE COMPONENT FLUID WITH PHASE CHANGE

In this Appendix we sketch a similar development for the case of a single component fluid with phase change, i.e., a volatile droplet in its own vapor. The van der Waals equation of state would be a first example.

Instead of (16) and (17), we will need to begin from the full compressible equations, including the energy equation

$$ \frac{\partial \rho}{\partial t} + (\rho u_i), \quad i = 0, \qquad (C1) $$

$$ \frac{\partial \rho u_i}{\partial t} + (\rho u_i u_j),_{j} = -P_{ij,j} + (\mu(u_{i,j} + u_{j,i})),_{j}, \quad (C2) $$

$$ \frac{\partial \rho e}{\partial t} + (\rho e u_j),_{j} = -P_{jk}u_{k,j} + \mu(u_{k,j} + u_{j,k})u_{k,j} + (\lambda T_j),_{j}, \quad (C3) $$

$$ P_{ij} = \left( p(\rho, T) - \kappa \left( \rho \rho_{,kk} + \frac{1}{2} \rho_{,k} \rho_{,k} - \frac{\rho}{T} \rho_{,k} T_{,k} \right) \right) \delta_{ij} + \kappa \rho_{,i} \rho_{,j}. \quad (C4) $$

In addition to density, velocity, and viscosity $ρ$, $u_i$, and $μ$, we have internal energy $e$, temperature $T$, and heat conductivity $λ$, $μ$ and $λ$ taken constant for simplicity. $P_{ij}$ denotes the pressure tensor, including the Korteweg stress [41], and $κ$ is the gradient energy coefficient. The function $p(ρ, T)$ is the equation of state.

The equilibrium solution to these equations, corresponding to (35), is found for $u_i = 0$ and $T = \text{const}$, and $P_{ij} = p_s(T)\delta_{ij}$, where $p_s(T)$ is the liquid-vapor coexistence pressure at temperature $T$. In implicit form this can be written as

$$ f(\rho, T) + p_s(T) = \frac{\kappa}{2} \rho_x^2, \quad (C5) $$

where the free energy $f$ is related to the pressure as $p = -f + \rho f_ρ$ (indices denoting derivatives).

Analogous to (31)–(38), we now examine the departure from this equilibrium due to interface curvature, and thus derive the Gibbs-Thomson relation for this specific diffuse interface model with a gradient energy. We make a change of variable as before, $r = R + \epsilon x$ where $\epsilon \ll 1$ indicates that curvature is much larger than the interface width. We introduce the perturbed quantities as

$$ \begin{aligned} P_{ij} &= p_s(T)\delta_{ij} + \epsilon P_{ij}^1, \\ \rho &= \rho_0 + \epsilon \rho_1, \\ u_i &= 0 + \epsilon u_i^1, \\ e &= e_0 + \epsilon e_1, \\ T &= T_s + \epsilon T_1. \end{aligned} \quad (C6) $$

Introducing these in $P_{11}$ and expanding we find

$$ \begin{aligned} P_{11} &= P_{11}^0 + \epsilon P_{11}^1, \\ P_{11}^0 &= p_s(T_s), \\ P_{11}^1 &= T_1 \frac{d p_s}{dT} + \kappa (\rho_1 \rho_0 f_{\rho\rho} + \rho_{0x} \rho_{1x} - \rho_1 \rho_{0xx} - \rho_0 \rho_{1xx}) \\ &\quad + \kappa \frac{1}{R} \rho_0 \rho_{0x}. \end{aligned} \quad (C7) $$

Multiplying $P_{11}^1$ by $\rho_{0x}/\rho_0^2$ and integrating we can obtain a solvability condition similar to (38). The second term in $P_{11}^1$ can be shown to vanish using the relation in (C5). The remaining terms are then

$$ T_1 \frac{d p_s}{dT} \int_v^l \frac{\rho_{0x}}{\rho_0^2} dx + \kappa \frac{1}{R} \int_v^l \frac{\rho_{0x}^2}{\rho_0} dx = 0, $$
where

$$ \int_v^l \frac{\rho_{0x}}{\rho_0^2} dx = \frac{1}{\rho_v} - \frac{1}{\rho_l} = \frac{1}{\tilde{\rho}}, \\ \kappa \int_v^l \frac{\rho_{0x}^2}{\rho_0} dx = \int_v^l \frac{\sqrt{2\kappa(f+p_s)}}{\rho_0} d\rho = \frac{\sigma}{\tilde{\rho}}. \quad (C8) $$

Here we denote the density corresponding to difference in specific volume $\tilde{\rho}$ as $\tilde{\rho} = 1/(1/\rho_v - 1/\rho_l)$. In the last integral the equilibrium solution in (C5) was used, we identified the surface tension as $\sigma = \int_v^l \sqrt{2\kappa(f+p_s)} d\rho$, and introduced an averaged density as

$$ \bar{\rho} = \frac{\int_v^l \sqrt{2\kappa(f+p_s)} d\rho}{\int_v^l \frac{\sqrt{2\kappa(f+p_s)}}{\rho_0} d\rho}. \quad (C9) $$

The Clapeyron equation in this context reads

$$ \frac{dp_s}{dT} = \frac{L}{T(1/\rho_v - 1/\rho_l)} = \frac{L\tilde{\rho}}{T}. \quad (C10) $$

Here L denotes the latent heat. The first equation in (C8) can then be written as

$$ T_1 \frac{L}{\bar{T}} + \frac{1}{R} \frac{\sigma}{\bar{\rho}} = 0, \quad (C11) $$

i.e., the Gibbs-Thomson relation. We note the explicit expression for the average density in (C9), and that this would emphasize the low density of the gas phase.

To proceed, we now derive a scaling for the region in the vicinity of the contact line, with a length scale $\delta$, to be determined. The velocity scale $u$ is the liquid velocity. The lowest order temperature is the equilibrium temperature $T_s$, and $T_1$ is the perturbation of the temperature. In analogy with (39) we assume that the main balance in the energy equation (C3) is between the thermal diffusion $(\lambda T, j)_j$ and the $pdV$ like term $P_{jk}u_{k,j}$, which at first order is $p_s(\rho_l/\rho_v)(u/\delta)$.

Depending on the fluid properties it is also possible that the thermal diffusion is balanced instead by the convective term on the left hand side of (C3), which would give slightly different estimates. As in (39) we introduce a nondimensional order unity factor $b$, which is to be calibrated by comparison between the model and full simulations

$$ p_s \frac{\rho_l}{\rho_v} \frac{bu}{\delta} = \lambda \frac{T_1}{\delta^2}. \quad (C12) $$

In the momentum equation (C2) we take the balance between the viscous term and the surface tension, estimating $P_{11}$ as $P_{11} \sim T_1(dp_s/dT)$ and using (C10). As in (40) we introduce a nondimensional order unity parameter $a$

$$ T_1 \frac{L\tilde{\rho}}{T_s} \frac{1}{\delta} = \mu \frac{au}{\delta^2}. \quad (C13) $$

Here $1/\tilde{\rho} = (1/\rho_v - 1/\rho_l)$.

Equations (C11)–(C13) are now the counterparts to Eqs. (38)–(40), and are to be solved for $\delta$, $R$ and $T_1$. The result of this, corresponding to Eqs. (41), (42), and (38), are

$$ R = -\frac{1}{a} \frac{\delta}{(\mu u / \sigma)} \frac{\tilde{\rho}}{\bar{\rho}} = -\frac{1}{\sqrt{ab}} \frac{1}{(\mu u / \sigma)} \frac{\tilde{\rho}}{\bar{\rho}} \sqrt{\frac{\rho_v}{\rho_l}} \sqrt{\frac{\lambda \mu}{\tilde{\rho} p_s (L/T_s)}}, \quad (C14) $$

$$ \delta = \sqrt{\frac{a}{b}} \sqrt{\frac{\rho_v}{\rho_l}} \sqrt{\frac{\lambda \mu}{p_s \tilde{\rho} (L/T_s)}}, \quad (C15) $$

$$ T_1 = u \sqrt{ab} \sqrt{\frac{p_s \mu}{\lambda(L/T_s)} \frac{\rho_l}{\rho_v \tilde{\rho}}}. \quad (C16) $$

We can now write down the model for the apparent contact angle, as in (43), which is seen to be very similar:

$$ \cos(\theta_a) = \cos(\theta_e) + \frac{\delta}{R} = \cos(\theta_e) + a \frac{\tilde{\rho}}{\bar{\rho}} \text{Ca}, \quad (C17) $$

with the capillary number $Ca = \mu u / \sigma$.

[1] T. D. Blake and Y. D. Shikhmurzaev, Dynamic wetting by liquids of different viscosity, J. Colloid Interface Sci. 253, 196 (2002).

[2] L. M. Hocking, Rival contact-angle models and the spreading of drops, J. Fluid Mech. 239, 671 (1992).

[3] J. Eggers and R. Evans, Comment on “Dynamic wetting by liquids of different viscosity,” by T.D. Blake and Y.D. Shikhmurzaev, J. Colloid Interface Sci. 280, 537 (2004).

[4] M. Wilson, J. Summers, Y. D. Shikhmurzaev, A. Clarke, and T. D. Blake, Nonlocal hydrodynamic influence on the dynamic contact angle: Slip models versus experiment, Phys. Rev. E 73, 041606 (2006).

[5] Y. D. Shikhmurzaev, Singularities at the moving contact line: Mathematical, physical and computational aspects, Physica D 217, 121 (2006).

[6] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, Wetting and spreading, Rev. Mod. Phys. 81, 739 (2009).

[7] J. H. Snoeijer and B. Andreotti, Moving contact lines: Scales, regimes, and dynamical transitions, Annu. Rev. Fluid Mech. 45, 269 (2013).

[8] C. Huh and L. E. Scriven, Hydrodynamic model of steady movement of a solid/liquid/fluid contact line, J. Colloid Interface Sci. 35, 85 (1971).

[9] U. Lācis, P. Johansson, T. Fullana, B. Hess, G. Amberg, S. Bagheri, and S. Zaleski, Steady moving contact line of water over a no-slip substrate, Eur. Phys. J.: Spec. Top. 229, 1897 (2020).

[10] T. Fullana, S. Zaleski, and S. Popinet, Dynamic wetting failure in curtain coating by the volume-of-fluid method, Eur. Phys. J.: Spec. Top. 229, 1923 (2020).

[11] D. Gründing, M. Smuda, T. Antritter, M. Fricke, D. Rettenmaier, F. Kummer, P. Stephan, H. Marschall, and D. Bothe, A comparative study of transient capillary rise using direct numerical simulations, Appl. Math. Model. 86, 142 (2020).

[12] S. Afkhami and M. Bussmann, Height functions for applying contact angles to 3D VOF simulations, Int. J. Numer. Meth. Fluids 61, 827 (2009).

[13] D. Legendre and M. Maglio, Comparison between numerical models for the simulation of moving contact lines, Comput. Fluids 113, 2 (2015).
[14] P. D. M. Spelt, A level-set approach for simulations of flows with multiple moving contact lines with hysteresis, *J. Comput. Phys.* **207**, 389 (2005).

[15] Y. Kulkarni, T. Fullana, and S. Zaleski, Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized Navier boundary condition, *Proc. R. Soc. A* **479**, 20230141 (2023).

[16] T. Fullana, Y. Kulkarni, M. Fricke, S. Popinet, S. Afkhami, D. Bothe, and S. Zaleski, A consistent treatment of dynamic contact angles in the sharp-interface framework with the generalized Navier boundary condition, *arXiv:2411.10762*.

[17] D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, *J. Comput. Phys.* **155**, 96 (1999).

[18] P. Yue, C. Zhou, and J. J. Feng, Sharp-interface limit of the Cahn-Hilliard model for moving contact lines, *J. Fluid Mech.* **645**, 279 (2010).

[19] U. Lācis, M. Pellegrino, J. Sundin, G. Amberg, S. Zaleski, B. Hess, and S. Bagheri, Nanoscale sheared droplet: Volume-of-fluid, phase-field and no-slip molecular dynamics, *J. Fluid Mech.* **940**, A10 (2022).

[20] R. Scardovelli and S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, *Annu. Rev. Fluid Mech.* **31**, 567 (1999).

[21] S. Popinet and S. Zaleski, A front-tracking algorithm for accurate representation of surface tension, *Int. J. Numer. Methods Fluids* **30**, 775 (1999).

[22] G. Tryggvason, R. Scardovelli, and S. Zaleski, *Direct Numerical Simulations of Gas-Liquid Multiphase Flows* (Cambridge University Press, Cambridge, 2011).

[23] https://basilisk.fr/.

[24] S. Popinet, Gerris: A tree-based adaptive solver for the incompressible euler equations in complex geometries, *J. Comput. Phys.* **190**, 572 (2003).

[25] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, *J. Comput. Phys.* **228**, 5838 (2009).

[26] S. Popinet, A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations, *J. Comput. Phys.* **302**, 336 (2015).

[27] J. U. Brackbill, D. B. Kothe, and C. Zemach, A continuum method for modeling surface tension, *J. Comput. Phys.* **100**, 335 (1992).

[28] J. B. Bell, P. Colella, and H. M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, *J. Comput. Phys.* **85**, 257 (1989).

[29] S. Afkhami and M. Bussmann, Height functions for applying contact angles to 2D VOF simulations, *Int. J. Numer. Meth. Fluids* **57**, 453 (2008).

[30] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, *J. Chem. Phys.* **28**, 258 (1958).

[31] D. Jacqmin, Contact-line dynamics of a diffuse fluid interface, *J. Fluid Mech.* **402**, 57 (2000).

[32] A. Carlson, Capillarity and dynamic wetting, Ph.D. thesis, KTH Royal Institute of Technology, 2012.

[33] T.-Z. Qian, X.-P. Wang, and P. Sheng, Molecular scale contact line hydrodynamics of immiscible flows, *Phys. Rev. E* **68**, 016306 (2003).

[34] F. Hecht, New development in freefem++, *J. Numer. Math.* **20**, 251 (2012).

[35] O. Pironneau, F. Hecht, A. L. Hyaric, and J. Morice, Freefem++ (2020), https://freefem.org/.

[36] https://github.com/UgisL/FreeFEM-NS-CH.

[37] T.-Z. Qian, X.-P. Wang, and P. Sheng, A variational approach to moving contact line hydrodynamics, *J. Fluid Mech.* **564**, 333 (2006).

[38] A. J. Briant and J. M. Yeomans, Lattice Boltzmann simulations of contact line motion. II. Binary fluids, *Phys. Rev. E* **69**, 031603 (2004).

[39] T. Fullana, Simulation and optimization of complex phenomena in multiphase flows, Ph.D. thesis, Sorbonne University, 2022.

[40] T. Fullana, S. Zaleski, and G. Amberg, Dataset for “Mass diffusion and bending in dynamic wetting by phase-field and sharp interface models” (2025), 10.5281/zenodo.16927267.

[41] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interface methods in fluid mechanics, *Annu. Rev. Fluid Mech.* **30**, 139 (1998).