empty

# Contact Angle Hysteresis on Rough Surfaces Part I: Mechanical Energy Balance
Framework

Dalton J.E. Harvie1<EMAIL_ADDRESS>1Department of Chemical
Engineering, University of Melbourne, Parkville, VIC, 3010, Australia

###### Abstract

Using as a starting point conservation of momentum, a multiphase mechanical
energy balance equation is derived that accounts for multiple material phases
and interfaces present within a moving control volume. This balance is applied
to a control volume that is anchored to a three phase contact line as it
advances over the surface of a rough and chemically homogeneous solid. Using
semi-quantitative models for the material behaviour occurring within the
control volume, an order-of-magnitude analysis is performed to find what terms
within the balance are significant, producing an equation that can be used to
predict contact angle hysteresis from a knowledge of interface dynamics
occurring around the three phase contact line. In addition to this equation,
the theory also answers several questions that have been discussed within the
wetting literature: Namely that (static) contact angle hysteresis is a
function of conditions around the three phase contact line, as opposed to the
surrounding flow system; That contact angle hysteresis results from interface
‘jumps’ that dissipate energy, rather than directly from contact line
deformation; That interfacial dynamics is required to predict contact angle
hysteresis, but that these dynamics should be interpreted via energy
conservation, and; That dynamic contact angles depend on kinetic energy
transport around the three phase contact line, as well as local energy
dissipation. The framework has been derived using assumptions of
incompressible Newtonian fluids, reversible interface formation and zero-
strain solids — future work could relax these assumptions to make the theory
more generally applicable.

###### keywords:

wetting, contact angle hysteresis, mechanical energy, rough surfaces

##  1 Introduction

The ability of a liquid to ‘wet’ a solid is described by the angle that the
liquid makes with the solid when the interface is stationary — the static
contact angle. For a liquid in contact with an ideal, smooth and chemically
homogeneous surface, this angle $\theta_{\mathrm{e}}$ is unique and given by
Young’s Equation. Real surfaces, which are often rough on a variety of length
scales, are chemically heterogeneous and/or involve some type of irreversible
work of adhesion, display a range of equilibrium contact angles: The maximum
is the static advancing angle $\theta_{\mathrm{a}}$, above will the interface
will advance, and the minimum the static receding angle,
$\theta_{\mathrm{r}}$, below which it will recede. The difference between
these two angles is defined as the range of contact angle hysteresis (CAH).
CAH, $\theta_{\mathrm{a}}$ and $\theta_{\mathrm{r}}$ are critical wetting
parameters that determine (for example) how easily drops can move over solid
surfaces, under what conditions liquid films will smoothly coat surfaces or
whether gas injection will aid particle floatation. Technologies that depend
on CAH angles include established processes such as industrial separation
devices or the wetting behaviour of fabrics, through to more novel processes
such as transparent self-cleaning surfaces for solar power generation or low-
energy liquid fuel separation membranes (Cassie & Baxter, 1944; Wu et al.,
2002; Feng et al., 2004; Sun et al., 2005; Callies & Quéré, 2005; Li et al.,
2013). For the design and optimisation of these processes general and
validated wetting theories are needed that can predict the CAH range. However,
as highlighted via several recent works such theories are not yet available,
with fundamental questions remaining about the nature of the wetting
process(Eral et al., 2013; Jiang et al., 2019; Butt et al., 2022). The purpose
of this study is to derive an energy conservation framework that can be
applied to predict CAH angles.

The early energy-based wetting theories of Wenzel and Cassie remain
influential in interpreting wetting phenomena. Considering the energy change
that occurs as a liquid/gas interface advances a small distance over the
surface of a rough solid, Wenzel (1936) proposed that the apparent contact
angle is related to the roughness $r$ of the surface, defined as the total to
projected surface area ratio (see section 4.1). In deriving this theory,
Wenzel assumed that the liquid completely wets each surface undulation. Cassie
& Baxter (1944) recognised that such ‘complete’ wetting did not necessarily
occur, and derived an expression for the apparent contact angle on a partially
wet surface in terms of the wetted and non-wetted liquid areas per projected
solid area, being $f_{1}$ and $f_{2}$, respectively (see section 4.1). The
Wenzel and Cassie & Baxter equations are useful for understanding and
interpreting experimental data, however, theories based solely on these
concepts (e.g. Bico et al., 2002; Patankar, 2003) are not predictive as the
proportion of solid surface wetted by a liquid is not known _a priori_. Also,
these theories give only one static contact angle for a rough structured
surface, rather than the CAH range that is observed experimentally.

Another series of works is based on the concept of energy minimisation of an
entire drop sitting on a rough surface. Johnson Jr & Dettre (1964) computed
the free energy of drops residing at the centre of concentric sinusoidal
roughness rings, showing that the energy of the system oscillated as the drop
volume increased and the interface advanced over each ring. They interpreted
the amplitude of these energy oscillations as energy barriers that must be
overcome by macroscopic vibrational energy to allow interface movement,
implying that as the height of surface roughness decreases, the range of CAH
should also decrease. This conclusion is contrary to experimental evidence
however which shows that roughness-induced CAH strongly depends on surface
topology (relative shape), rather than absolute roughness size, provided that
gravitational and Laplace pressure effects can be neglected on the lengthscale
of the roughness (Öner & McCarthy, 2000; Dorrer & Rühe, 2008; Li et al., 2016;
Jiang et al., 2019). Other studies have used similar static free-energy
minimisation concepts to explain CAH for a variety of periodically shaped
surfaces(Extrand, 2002; Brandon et al., 2003; Marmur, 2006, 2022), however in
general the results do not qualitatively agree with observation. For example,
Brandon et al. (2003) used minimal surface energy modelling to show that the
apparent contact angle range for a drop on a doubly periodic undulating
surface approached a single value (the Cassie angle) as the drop size to
roughness ratio increased, again contradicting the above referenced
experimental observations that show that the CAH range becomes quite constant
at large droplet to roughness size ratios. A related question also remains
about these analyses: Is it really necessary to consider the energy of the
entire flow system (most commonly a droplet) to calculate the CAH range, or is
it instead a property associated with the three phase contact line (TPCL) that
can be applied to a wide variety of surrounding flow systems? This is an
unresolved question that has garnered conflicting opinions (Gao & McCarthy,
2007; McHale, 2007; Gao & McCarthy, 2007; Nosonovsky, 2007; Panchagnula &
Vedantam, 2007; Marmur, 2022).

Other studies conclude that CAH is substantially a property of the TPCL
region(Nosonovsky, 2007; Panchagnula & Vedantam, 2007) and that CAH has its
origin in the energy dissipation that occurs around the TPCL as it advances
over the rough solid. Central to this concept is that the advancing and
receding angles are defined as those measured while the interface is moving
(albeit at a vanishingly slow velocity) rather than being determined solely by
static thermodynamic states. An influential study in this vein is Joanny & de
Gennes (1984) who proposed a model for CAH on a surface that contains a dilute
number of ‘strong defects’ as a model for surface contamination or dilutely
distributed surface roughness. The theory considered the ‘pinning’ and
subsequent ‘depinning’ or ‘jumping’ of the contact line as it advanced over a
surface, assuming that during each interface jump surface potential energy is
dissipated to heat. Joanny & de Gennes calculated this energy dissipation
amount under ideal conditions and incorporated it into an equation for CAH
relevant to dilute defect surfaces. More recent works have experimentally
observed the pinning/depinning behaviour of the fluid interface near the
TPCL(Jiang et al., 2019; Priest et al., 2009; Forsberg et al., 2010). Other
studies have used a variety of energy conservation principles to extend Joanny
& de Gennes’s work to periodic surfaces(Raj et al., 2012; Butt et al., 2017;
Jiang et al., 2019) or interpreted measurements of CAH in terms contact line
energy dissipation and interfacial ‘jumping’ dynamics(Priest et al., 2007,
2013; Dorrer & Rühe, 2008; Song et al., 2022). Despite these successes however
questions remain about this conceptual model of CAH; around what specific TPCL
region should energy be conserved, how should the energies of real rough
surfaces that may contain randomly shaped structures or micro bubbles/droplets
be incorporated in the energy analysis, and why must the energy dissipated due
to interface jumping dynamics be balanced by only interfacial (rather than
material) stress movements?

On a slightly different track Joanny & de Gennes (1984)’s work also explored
the influence of ’weak’ surface defects on CAH, referring to smooth defects as
those that cause the fluid interface near the solid to become distorted, but
that do not result in the aforementioned ‘pinning’ and ‘depinning’ behaviour
of the TPCL. Joanny & de Gennes’s conclusion was that isolated weak defects do
not generally result in hysteresis, however other works have extended this
analysis to conclude that distributions of weak surface heterogeneities can
cause CAH(Pomeau & Vannimenus, 1985; Robbins & Joanny, 1987; Öpik, 2000). In
related work Cox (1983) examined how an interface changes as it moves over a
gently undulating sinusoidal periodic rough surface, showing that when the
interface moved in the direction of roughness periodicity TPCL ‘jumps’
occurred (i.e., ‘strong’ defects leading to CAH), but when advancing in other
directions relative to the periodicity direction the interface moved
continuously (i.e., ‘weak’ defects producing no CAH). It should be noted that
these theoretical studies predict the possible shapes that a fluid interface
can take when passing over arrangements of surface heterogeneities, and from
these the range of CAH angles are inferred: In general the link between CAH
angles calculated via these interface topology methods and those calculated by
energy conservation has not been established.

Reviewing this body of literature key questions about wetting behaviour
remain. Specifically there is conjecture regarding how energy conservation,
fluid interface distortion and the dynamics of interface ‘jumping’ can be
combined to predict CAH. Critically, there is no experimental consensus or
fundamental analysis that shows how energy conservation should be applied to
predict CAH — works that are based on energy conservation around entire
droplets in general do not predict experimentally observed CAH trends, while
studies that are based on energy changes around a moving TPCL lack a rigorous
energy conservation basis. Further, while many studies do view contact line
jumping as a source of CAH (‘strong’ defects), there is confusion over whether
fluid interface distortion (‘weak’ defects) in isolation can produce CAH, and
how studies that predict CAH angles via interface topology and dynamics can be
mathematically related to CAH angles resulting from energy conservation.

This study addresses these questions. Specifically we derive a rigorous energy
conservation framework which when applied to the moving TPCL can be used to
predict contact angle hysteresis (CAH). We do this by first deriving a general
multiphase mechanical energy balance (MMEB) equation that (section 2) is
equivalent to the conventional single phase mechanical energy balance equation
but includes terms associated with interfacial stresses acting between each of
the material phases. We then apply this MMEB equation to a cylindrical control
volume that is anchored to a TPCL as it moves across a rough solid surface at
a vanishingly slow speed (section 3). By adopting order-or-magnitude models
that describe how each of the materials behave within the control volume
during the advance, we neglect terms that are insignificant to the calculation
and derive a resulting energy conservation equation that predicts CAH given
knowledge of the interfacial dynamics occurring within the TPCL. A strength of
the analysis is that mathematical constraints are derived that determine under
what physical conditions the theory will be valid (section 4), giving insight
into past and future modelling theories and experimental wetting studies.

##  2 Macroscopic Mechanical Energy Balance for Immiscible Multiphase Mixtures

In this section we derive a macroscopic mechanical energy balance for a moving
control volume that contains a number of immiscible phases (see figure 1). The
derivation is similar to that of Bird et al. (2002, p221, §7.8) except that
interfacial tension acts at the interface between each pair of phases, and the
balance is not specific to Newtonian liquids111We do assume that each material
has a symmetric stress tensor, however for most homogeneous materials this
assumption appears to be valid (Kuiken, 1995; Dahler & Scriven, 1961)..

$\boldsymbol{v}_{\mathrm{CV}}$$V_{\mathrm{CV}}$$\boldsymbol{n}_{\mathrm{S},23}$$\boldsymbol{n}_{\mathrm{S},23}$$\boldsymbol{n}_{\mathrm{S},13}$$\boldsymbol{n}_{\mathrm{CV}}$$\boldsymbol{n}_{\mathrm{S},12}$$\boldsymbol{n}_{\mathrm{S},34}$$S_{\mathrm{CV}}$

Figure 1: The macroscopic mechanical energy balance is performed on a moving
control volume that contains multiple immiscible phases bounded by multiple
interface types. Each interface type has an associated unit normal vector
$\boldsymbol{n}_{\mathrm{S},ij}$ that is directed into phase $i$ (where
$i<j$). The velocity and outward unit normal of the control volume boundary
are $\boldsymbol{v}_{\mathrm{cv}}$ and $\boldsymbol{n}_{\mathrm{cv}}$,
respectively. In the indicated configuration four material phases are present
within the control volume, and they intersect along four different interface
types.

Our starting point is a momentum equation for a mixture of immiscible phases,

$\frac{\partial}{\partial
t}\rho\boldsymbol{v}_{\mathrm{}}+\boldsymbol{\nabla}_{\mathrm{}}\cdot\rho\boldsymbol{v}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}=\boldsymbol{\nabla}_{\mathrm{}}\cdot[\mathsf{T}_{\mathrm{M}}+\mathsf{T}_{\mathrm{S}}]+\rho\boldsymbol{g}_{\mathrm{}}.$
(1)

Here $\boldsymbol{v}_{\mathrm{}}$ is the local phase velocity (assumed to vary
continuously throughout), $\rho$ the local phase density,
$\mathsf{T}_{\mathrm{M}}$ the local total material stress at any point within
any phase and $\mathsf{T}_{\mathrm{S}}$ the local surface stress acting on the
interfaces between phases. The functional form of the material stress tensor
$\mathsf{T}_{\mathrm{M}}(\boldsymbol{x}_{\mathrm{}})$ is a property of the
material type present at $\boldsymbol{x}_{\mathrm{}}$.

Formally, for equation (1) to be valid everywhere within $V_{\mathrm{cv}}$ all
terms appearing in the equation must be defined not only within each phase but
also on the interfaces between phases. This includes terms such as
$\mathsf{T}_{\mathrm{M}}$ and $\rho$ that are associated with a particular
material type. Such formal definitions could be made; however as equation (1)
is integrated over space in the following analysis, as long as any phase-
specific terms are finite on each interface then their interface values do not
affect the final energy balance. Hence, we simply assume that the interface
values for $\mathsf{T}_{\mathrm{M}}(\boldsymbol{x}_{\mathrm{}})$ and $\rho$
are finite.

The effects of interfacial tension on material movement are captured using a
surface stress tensor, $\mathsf{T}_{\mathrm{S}}(\boldsymbol{x}_{\mathrm{}})$.
The form of this tensor is taken from Lafaurie et al. (1994), but generalised
here to include interfaces between multiple phases;

$\mathsf{T}_{\mathrm{S}}(\boldsymbol{x}_{\mathrm{}})=\sum_{i<j}\sigma_{ij}(\mathsf{I}_{\mathrm{}}-\boldsymbol{n}_{\mathrm{S},ij}\boldsymbol{n}_{\mathrm{S},ij})\delta_{\mathrm{S},ij}.$
(2)

In this equation $\sigma_{ij}$ is the constant surface energy per unit area
(or interfacial tension) associated with the ‘$ij$ interface’ (i.e., the
interface between phases $i$ and $j$), $\mathsf{I}_{\mathrm{}}$ is the
identity tensor, $\boldsymbol{n}_{\mathrm{S},ij}$ is a unit vector directed
normal to the $ij$ interface and into phase $i$, and $\delta_{\mathrm{S},ij}$
is a ‘surface’ delta function that is nonzero only on the $ij$ interface. The
surface delta function is essentially a multidimensional analogue of the Dirac
delta function and has been utilised extensively in the development of
computational fluid dynamics methods (Brackbill et al., 1992; Lafaurie et al.,
1994). In the present context it has the property that

$\int_{V_{\mathrm{cv}}}\delta_{\mathrm{S},ij}dV=A_{ij},$ (3)

where $A_{ij}$ is the total area of the $ij$ interface existing within the
control volume $V_{\mathrm{cv}}$. Further, under conditions where
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{n}_{\mathrm{S},ij}$ is uniform
over a particular intersection between a surface $S_{\mathrm{cv}}$ and an
interface defined by $\delta_{\mathrm{S},ij}$, the surface integral of the
surface delta function is given by

$\int_{S_{\mathrm{cv}}}\delta_{\mathrm{S},ij}dS=\frac{l_{ij}}{\sqrt{1-(\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{n}_{\mathrm{S},ij})^{2}}},$
(4)

where $l_{ij}$ represents the line length of the intersection between the two
surfaces and $\boldsymbol{n}_{\mathrm{cv}}$ is a unit normal to the surface
$S_{\mathrm{cv}}$. These properties of the surface delta function and others
are discussed in more detail in the Appendix A.

Note that in general interfacial stresses will act at each one of the
immiscible material boundaries that exist within $V_{\mathrm{cv}}$.
Consequently, the sum in equation (2) cycles through all possible phase
combinations under the condition that $i<j$; that is, $j=1\rightarrow m$ and
$i=1\rightarrow j$ where $m$ is the total number of material phases present.
Thus, stresses from a possible $(m-1)!$ interface types may be included in the
momentum balance.

By using equation (2) to represent surface stresses, three assumptions about
the system are implied. Firstly, as the surface stress is a sum of
contributions from each interface type, we have neglected any ‘line tension’
stresses that may occur at the intersection between interfaces. While no
concensus regarding the existence these stresses has been reached in the
literature, most studies suggest that even if line tension does exist, it has
a negligible effect on macroscopically measureable contact angles (Boruvka &
Neumann, 1977; Marmur, 1997; Pompe et al., 1999; Marmur, 2006). Secondly, by
assuming constant surface energies for each interface type, we have neglected
any Marangoni forces that would exist if surfactants or thermal gradients were
present within the control volume. Thirdly, by representing the surface stress
by equation (2) we have implicitly assumed that the process of surface
creation or destruction is reversible on a molecular scale. We discuss
implications of this assumption in section 4.3.

With the immiscible multiphase momentum equation defined, we proceed by taking
the dot product of equation (1) with the local velocity
$\boldsymbol{v}_{\mathrm{}}$ and then integrating the result over the volume
$V_{\mathrm{cv}}$. Noting that both the stress tensors
$\mathsf{T}_{\mathrm{M}}$ and $\mathsf{T}_{\mathrm{S}}$ are symmetric,
application of the Leibnitz formula for differentiating a volume integral,
Gauss-Ostrogradskii theorem and compressible continuity equation yields,

$\frac{d}{dt}\int_{V_{\mathrm{cv}}}\left(\frac{1}{2}\rho
v^{2}+\rho\hat{\Phi}\right)dV=\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\left(\frac{1}{2}\rho
v^{2}+\rho\hat{\Phi}\right)\left(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}}\right)\right]dS\\\
+\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\mathsf{T}_{\mathrm{M}}\cdot\boldsymbol{v}_{\mathrm{}}\right]dS-\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV\\\
+\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\mathsf{T}_{\mathrm{S}}\cdot\boldsymbol{v}_{\mathrm{}}\right]dS-\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{S}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV.$
(5)

Here $\boldsymbol{v}_{\mathrm{cv}}$ and $\boldsymbol{n}_{\mathrm{cv}}$ are the
velocity and outwardly directed unit normal of the control volume boundary
$S_{\mathrm{cv}}$, respectively, $v$ is the magnitude of the local velocity
$\boldsymbol{v}_{\mathrm{}}$, and $\hat{\Phi}$ is a conservative gravitational
potential function satisfying
$\boldsymbol{g}_{\mathrm{}}=-\boldsymbol{\nabla}_{\mathrm{}}\hat{\Phi}$.

To simplify equation (5) further we concentrate on the last two terms on the
right hand side which relate to interfacial stresses. For the first of these
we substitute in the surface stress definition of equation (2) to find

$\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\mathsf{T}_{\mathrm{S}}\cdot\boldsymbol{v}_{\mathrm{}}\right]dS=\sum_{i<j}\sigma_{ij}\int_{S_{\mathrm{cv}}}\delta_{\mathrm{S},ij}\left(\mathsf{I}_{\mathrm{}}-\boldsymbol{n}_{\mathrm{S},ij}\boldsymbol{n}_{\mathrm{S},ij}\right):\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dS$
(6)

For the second term we use the surface delta function transport equation
derived in the Appendix (section A),

$\frac{\partial\delta_{\mathrm{S},ij}}{\partial
t}+\boldsymbol{\nabla}_{\mathrm{}}\cdot\left(\delta_{\mathrm{S},ij}\boldsymbol{v}_{\mathrm{}}\right)=\delta_{\mathrm{S},ij}\left(\mathsf{I}_{\mathrm{}}-\boldsymbol{n}_{\mathrm{S},ij}\boldsymbol{n}_{\mathrm{S},ij}\right):\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}.$
(7)

Substituting $\mathsf{T}_{\mathrm{S}}$ from equation (2) into the second
interfacial stress term of equation (5), and then using the right hand side of
the equation (7) to expand the double dot product gives

$\displaystyle\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{S}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV$
$\displaystyle=\sum_{i<j}\int_{V_{\mathrm{cv}}}\sigma_{ij}\delta_{\mathrm{S},ij}\left(\mathsf{I}_{\mathrm{}}-\boldsymbol{n}_{\mathrm{S},ij}\boldsymbol{n}_{\mathrm{S},ij}\right):\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV$
$\displaystyle=\sum_{i<j}\int_{V_{\mathrm{cv}}}\frac{\partial(\sigma_{ij}\delta_{\mathrm{S},ij})}{\partial
t}+\boldsymbol{\nabla}_{\mathrm{}}\cdot\left(\sigma_{ij}\delta_{\mathrm{S},ij}\boldsymbol{v}_{\mathrm{}}\right)dV.$
(8)

Using the scalar Liebnitz theorem on the first term on the right of this
equation and the Gauss-Ostrogradskii theorem on the second term yields

$\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{S}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV=\\\
\sum_{i<j}\left\\{\frac{d}{dt}\left(\sigma_{ij}A_{ij}\right)-\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\sigma_{ij}\delta_{\mathrm{S},ij}(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}})\right]dS\right\\}$
(9)

where equation (3) has been used to relate surface area to the volume integral
of $\delta_{\mathrm{S},ij}$.

Finally, substituting equations (6) and (9) back into equation (5) and
simplfying the material stress surface integral gives the immiscible
multiphase mechanical energy balance valid for compressible and incompressible
materials,

$\underbrace{\frac{d}{dt}\left[\int_{V_{\mathrm{cv}}}\left(\frac{1}{2}\rho
v^{2}+\rho\hat{\Phi}\right)dV+\sum_{i<j}\sigma_{ij}A_{ij}\right]}_{\text{\parbox{135.97101pt}{\centering\scriptsize
rate of change of kinetic, gravitational potential and interfacial surface
energy within $V_{\mathrm{cv}}$\@add@centering}}}=\\\
\underbrace{\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\left[\left(\frac{1}{2}\rho
v^{2}+\rho\hat{\Phi}+\sum_{i<j}\sigma_{ij}\delta_{\mathrm{S},ij}\right)\left(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}}\right)\right]dS}_{\text{\parbox{168.27927pt}{\centering\scriptsize
rate at which kinetic, graviational potential and interfacial surface energy
are advected into $V_{\mathrm{cv}}$\@add@centering}}}\\\
+\underbrace{\sum_{i<j}\sigma_{ij}\int_{S_{\mathrm{cv}}}\delta_{\mathrm{S},ij}\left(\mathsf{I}_{\mathrm{}}-\boldsymbol{n}_{\mathrm{S},ij}\boldsymbol{n}_{\mathrm{S},ij}\right):\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dS}_{\text{\parbox{130.95467pt}{\centering\scriptsize
rate of work done on the contents of $V_{\mathrm{cv}}$ by interfacial tension
acting at $S_{\mathrm{cv}}$\@add@centering}}}\\\
+\underbrace{\int_{S_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dS}_{\text{\parbox{54.7338pt}{\centering\scriptsize
rate of work done on the contents of $V_{\mathrm{cv}}$ by material stresses
acting at
$S_{\mathrm{cv}}$\@add@centering}}}-\underbrace{\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}dV}_{\text{\parbox{53.72205pt}{\centering\scriptsize
rate at which energy dissipates to heat via material stresses acting within
$V_{\mathrm{cv}}$\@add@centering}}}$ (10)

Along with the usual terms found in the single phase mechanical energy balance
(Bird et al., 2002, p81, §3.3) equation (10) contains three interfacial stress
terms: The first represents the rate of change of interfacial energy contained
within the control volume; the second the rate at which interfacial energy is
advected across the control volume surface; and the third the rate at which
interfacial stresses perform work on the control volume at the control volume
surface. Within the next section we demonstrate how these terms are evaluated
for a specific control volume geometry.

##  3 Calculating the advancing contact angle on a rough solid surface

In this section we analyse the macroscopic contact angle of an interface
between two immiscible fluids that slowly advances over a rough solid. The
analysis uses the multiphase mechanical energy balance derived in section 2,
applied to a small control volume (CV) which moves with the advancing three
phase contact line across the surface of the solid. Using semi-quantitative
models for material and interface behaviour, an order of magnitude analysis is
performed to determine which terms within the energy balance are significant,
and from this an expression for the advancing contact angle is found that is
valid in the limit of an infinitely slowly moving interface.

###  3.1 Defining the physical system and moving control volume geometry

Figure 2 illustrates the physical system that is analysed. Two fluid phases,
labelled $\bigcirc$1 and $\bigcirc$2, are bounded below by a solid phase,
labelled $\bigcirc$S. As the two fluids are immiscible they are separated by a
distinct fluid interface. The fluids may be either liquids or gases, but for
this study are assumed to be completely insoluble with each other and with the
solid, incompressible, and above a certain lengthscale $h_{\mathrm{mol}}$ (for
molecular), act as continua. The implications of these assumptions are
discussed further below and in section 4. The surface of the solid is rough,
having undulations of a characteristic size $h_{\mathrm{rough}}$. A
hydrodynamic flow is occuring on a lengthscale of $h_{\mathrm{surround}}$
which is much larger than $h_{\mathrm{rough}}$. This flow slowly drives the
fluid interface to the right: Hence, phase $\bigcirc$1 is slowly advancing
over the solid while phase $\bigcirc$2 is slowly receding.

$h_{\mathrm{surround}}$$h_{\mathrm{rough}}$$r_{\mathrm{cv}}$controlvolume$\theta_{\mathrm{a}}$

(a)

control volumeboundaryaverage solidsurface planeinterface planeprojected
fluid$\underbrace{\rule[-4.30554pt]{0.0pt}{0.0pt}\,\sigma_{\mathrm{12}}\overleftarrow{A}_{\mathrm{12}}+\sigma_{\mathrm{1s}}\overleftarrow{A}_{\mathrm{1s}}+\sigma_{\mathrm{2s}}\overleftarrow{A}_{\mathrm{2s}}\,}_{\mbox{\rule[0.0pt]{0.0pt}{7.53471pt}$\overleftarrow{\sigma}A_{\mathrm{cv}}$}}$$\underbrace{\rule[-4.30554pt]{0.0pt}{0.0pt}\,\sigma_{\mathrm{1s}}\overrightarrow{A}_{\mathrm{1s}}+\sigma_{\mathrm{2s}}\overrightarrow{A}_{\mathrm{2s}}+\sigma_{\mathrm{12}}\overrightarrow{A}_{\mathrm{12}}\,}_{\mbox{\rule[0.0pt]{0.0pt}{7.53471pt}$\overrightarrow{\sigma}A_{\mathrm{cv}}$}}$$\boldsymbol{e}_{1}$$\sigma_{\mathrm{12}}$$\theta_{\mathrm{a}}$$\boldsymbol{v}_{\mathrm{cv}}$$\boldsymbol{e}_{3}$$D$

(b)

Figure 2: Frame (a) shows the region surrounding the control volume and
contact line on a macroscopic scale, while frame (b) shows the same region on
the scale of the control volume and solid surface roughness. The symbols
$\bigcirc$1, $\bigcirc$2, and $\bigcirc$S, indicate regions of phase $1$
fluid, phase $2$ fluid, and solid, respectively.

The multiphase mechanical energy balance is applied to a moving CV as it
advances at a constant velocity
$\boldsymbol{v}_{\mathrm{cv}}=v_{\mathrm{cv}}\boldsymbol{e}_{1}$ over a
distance of $X_{\mathrm{cv}}$ along the solid surface. This is illustrated in
figure 3. Note that $v_{\mathrm{cv}}$ is characteristic of the surrounding
flow. The CV contains and is approximately centred on the ‘three phase contact
line’ (TPCL), defined as the intersection between the advancing fluid
interface and rough solid surface. The moving CV has the geometry of a
cylinder with radius $r_{\mathrm{cv}}$ and length $l_{\mathrm{cv}}$. The
dimensions of the volume are smaller than that of the surrounding hydrodynamic
flow ($h_{\mathrm{surround}}$), yet larger than that of the solid surface
roughness ($h_{\mathrm{rough}}$). Hence, noting the above description of the
physical system and CV geometry we effectively assume the separation of four
lengthscales in our analysis,

$h_{\mathrm{mol}}\ll h_{\mathrm{rough}}\ll
r_{\mathrm{cv}},l_{\mathrm{cv}},X_{\mathrm{cv}}\ll h_{\mathrm{surround}}.$
(11)

Defining $\tau$ as the time taken for the CV to advance the entire distance
$X_{\mathrm{cv}}$, it follows that $X_{\mathrm{cv}}=v_{\mathrm{cv}}\tau$.

The precise centreline of the moving control volume is defined to lie at the
intersection between two averaged planes: the ‘average solid surface plane’
and ‘projected fluid interface plane’. The locations of these planes are
defined as those of the solid surface and fluid interfaces, respectively,
averaged over distances of $\mathcal{O}\left(r_{\mathrm{cv}}\right)$ (where
$\mathcal{O}\left(z\right)$ means ‘order $z$’). As the solid surface roughness
$h_{\mathrm{rough}}$ is of much smaller size than $r_{\mathrm{cv}}$, it
follows that the average solid surface plane is perfectly flat on the
lengthscale of the CV. In terms of the fluid interface, its topology is
governed by the momentum and surface stress equations (equations (1) and (2),
respectively), combined with boundary conditions specifying how the interface
interacts with the solid surface. The specific boundary condition that we
employ for the microscopic contact angle is Youngs equation, expressed as

$\left[\boldsymbol{n}_{\text{S},12}\right]_{\text{TPCL}}\cdot\boldsymbol{n}_{\mathrm{w}}=\cos\theta_{\mathrm{e}},$
(12)

where $\left[\boldsymbol{n}_{\text{S},12}\right]_{\text{TPCL}}$ is the unit
normal to the fluid interface at a point on the three phase contact line,
$\boldsymbol{n}_{\mathrm{w}}$ is the outwardly directed wall normal at the
same contact point and $\theta_{\mathrm{e}}$ is the equilibrium or ‘Youngs’
angle. As equation (12) specifies a direct relationship between the fluid
interface ($\boldsymbol{n}_{\text{S},12}$) and solid
($\boldsymbol{n}_{\mathrm{w}}$) normals along the TPCL, it follows that close
to the rough solid surface the fluid interface will have local curvatures that
are characteristic of the solid roughness — that is, of
$\mathcal{O}\left(1/h_{\mathrm{rough}}\right)$ — and that these curvatures
will exist within distances of $\mathcal{O}\left(h_{\mathrm{rough}}\right)$
from the TPCL. Conversely, further from the solid surface the topology of the
interface varies over the larger lengthscales of the surrounding flow (indeed,
this can be used to define $h_{\mathrm{surround}}$), so the curvature of the
interface there approaches $\mathcal{O}\left(1/h_{\mathrm{surround}}\right)$.
Hence, as $h_{\mathrm{rough}}\ll r_{\mathrm{cv}}\ll h_{\mathrm{surround}}$
(equation (11)), averaging the actual fluid interface over
$\mathcal{O}\left(r_{\mathrm{cv}}\right)$ produces a projected fluid interface
plane that is perfectly flat on the lengthscales of the CV, and as the
centreline of the control volume is defined as the intersection between the
projected fluid interface plane and average solid surface plane, on the scale
of the control volume its geometry is that of a perfect cylinder, with a
perfectly straight centreline.

This description of the fluid interface topology and geometry of the CV has
further implications for the mechanical energy balance application. At
distances of $\mathcal{O}\left(r_{\mathrm{cv}}\right)$ from the TPCL, the
actual fluid interface and projected fluid interface will at all times be
coincident and perfectly flat on these lengthscales. Hence, given that the
projected fluid interface is used to define the centreline of the CV, where
the actual fluid interface intersects the circumference of the CV (labeled as
$S_{\mathrm{cv,cir}}$ in Figure 3) it will be perfectly normal to the CV
boundary at all times. Similarly, along the ends of the CV (labelled as
$S_{\mathrm{cv,end}}$ in Figure 3) and at distances greater than
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$ from the TPCL, the actual fluid
interface will be flat, coincident with the projected fluid interface and
perfectly normal to the CV boundary. Closer to the TPCL however the interface
will undulate with curvatures of
$\mathcal{O}\left(1/h_{\mathrm{rough}}\right)$, crossing the CV boundary at
angles with are not necessarily normal to $S_{\mathrm{cv,end}}$.

averageplanesolid
surface$l_{\mathrm{cv}}$$x_{\mathrm{cv}}$$X_{\mathrm{cv}}$$A_{\mathrm{cv}}$$\boldsymbol{e}_{\mathrm{1}}$$\boldsymbol{e}_{\mathrm{2}}$$\boldsymbol{e}_{\mathrm{3}}$planefluid
interfaceprojected$S_{\mathrm{cv,end}}$$v_{\mathrm{cv}}$$S_{\mathrm{cv,cir}}$$\theta_{\mathrm{a}}$$S_{\mathrm{cv,end}}$$r_{\mathrm{cv}}$

Figure 3: The mechanical energy balance is performed over a cylindrical
control volume that is located at the intersection of the average solid
surface and projected fluid interface planes and moves forward over the solid
at a speed of $v_{\mathrm{cv}}$ and distance $X_{\mathrm{cv}}$.

The contact angle of phase $\bigcirc$1 on the macroscopic lengthscale is
defined as $\theta_{\mathrm{}}$. It is the angle between the average solid
surface plane and fluid interface plane, measured through phase $\bigcirc$1.
Equivalently, consistent with the above, it is the angle between the average
solid surface plane and actual fluid interface measured at distances of
$\mathcal{O}\left(r_{\mathrm{cv}}\right)$ from the TPCL. The objective of our
analysis is to determine the minimum macroscopic angle that _just_ causes the
fluid interface to advance continually over the solid surface, albeit at the
vanishingly slow velocity of $v_{\mathrm{cv}}$. An equivalent objective is to
find the maximum macroscopic angle that _just_ allows the fluid interface to
remain stationary. Either definition represents the advancing angle of phase
$\bigcirc$1 over solid $\bigcirc$S in the presence of phase $\bigcirc$2 and is
referred to as $\theta_{\mathrm{a}}$ (for ‘advancing’). Formally
$\theta_{\mathrm{a}}=\lim_{v_{\mathrm{cv}}\rightarrow 0}\theta$. Note that as
phase $\bigcirc$1 advances over the solid, phase $\bigcirc$2 recedes. Thus an
equivalent objective is to find the receding angle of phase $\bigcirc$2 (equal
to $\pi-\theta_{\mathrm{a}}$). Indeed, by swapping the physical properties
between phases $\bigcirc$1 and $\bigcirc$2 (detailed in section 4.1), we can
use the same analysis to determine the range of angles over which a fluid
interface will remain stationary — that is, the range of CAH.

###  3.2 Describing material dynamics within the moving control volume

As well as defining the physical system and CV geometry used in the mechanical
energy balance, to be able to perform an order of magnitude analysis on its
various terms we need to quantitatively describe how the materials within the
volume behave as a function of both space and time. Specifically we need
conceptual models for how the fluid velocities, pressures, interface topology
and solid stresses vary as the CV advances.

Within the fluid phase we assume that for the majority of the advancing time
$\tau$ the TPCL and surrounding fluid both move continuously at the slow
speeds of $\mathcal{O}\left(v_{\mathrm{cv}}\right)$. We refer to the system as
being in ‘equilibrium’ when this is the case and define the velocity field
existing during these times as
$\widetilde{\boldsymbol{v}_{\mathrm{}}}=\mathcal{O}\left(v_{\mathrm{cv}}\right)$.
However, at certain times during $\tau$, local areas of the TPCL will become
pinned by particular surface defects, creating local interface deflections
that become larger as the remainder of the TPCL continues to advance.
Eventually, once the surrounding TPCL has advanced some distance of
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$ from the pinning location, these
contact line sections will either detach from the surface defect and return to
the main TPCL, or detach from the TPCL and form isolated bubbles/droplets of
entrained fluid within the surface roughness. Either way, these detachment
processes cause the local TPCL and surrounding fluid to move at much faster
capillary driven speeds than the continuous CV advance speed
($v_{\mathrm{cv}}$), causing a viscous dissipation of energy. This slip-stick
dissipative motion has been described previously as the cause of contact angle
hysteresis(Joanny & de Gennes, 1984; Raj et al., 2012; Butt et al., 2017;
Jiang et al., 2019) and as discussed in the Introduction has also been
experimentally observed(Jiang et al., 2019; Priest et al., 2009; Forsberg et
al., 2010). In this study we define the local interface and fluid speeds
associated with these capillary driven events as
$\widehat{\boldsymbol{v}_{\mathrm{}}}=\mathcal{O}\left(v_{\mathrm{cap}}\right)$,
and the total time during which there is a dissipation event occurring within
the control volume as $\tau_{\mathrm{cap}}$. We further assume that the number
of defects within the control volume ($N$) is small enough and the capillary
velocity ($v_{\mathrm{cap}}$) large enough that only one dissipation event
occurs within the CV at any one time. With these assumptions, and for
convenience assuming that the analysis duration $\tau$ commences and finishes
while the system is in equilibrium, we can split the total time over which the
analysis is being conducted $\tau$ into a number of ‘dissipation events’ ($N$)
and ‘equilibrium stages’ ($N+1$), with the $k$th dissipation event starting at
$\widehat{t}_{k}$ and lasting for ${\widehat{\Delta t}}_{k}$, and the $k$th
equilibrium stage starting at $\widetilde{t}_{k}$ and lasting for
$\widetilde{\Delta t}_{k}$. The schematic timeline of Figure 4 illustrates
this decomposition. The following relationships result:

$\begin{gathered}\tau_{\mathrm{cap}}=\sum_{k=1}^{N}{\widehat{\Delta
t}}_{k}\quad\quad\tau=\sum_{k=1}^{N+1}\widetilde{\Delta
t}_{k}+\tau_{\mathrm{cap}}\\\ {\widehat{\Delta
t}}_{k}=\widetilde{t}_{k+1}-\widehat{t}_{k}\quad\quad\widetilde{\Delta
t}_{k}=\widehat{t}_{k}-\widetilde{t}_{k}\\\
\widetilde{t}_{1}=0\quad\quad\widetilde{t}_{N+1}+\widetilde{\Delta
t}_{N+1}=\tau\end{gathered}$ (13)

In terms of notation used in the remainder of the analysis, variables
annotated with a ‘hat’ correspond to properties associated with individual
dissipation events (where parts of the TPCL are moving at
$\mathcal{O}\left(v_{\mathrm{cap}}\right)$), variables annotated with a
‘tilde’ correspond to properties associated with the system while in
equilibrium (where the entire TPCL is moving at at most
$\mathcal{O}\left(v_{\mathrm{cv}}\right)$), and variables annotated with a
‘bar’ correspond to the entire advancing period $\tau$. The decomposition of
$\tau$ into separate dissipation and equilibrium stages is a key concept used
in the subsequent energy analysis.

$\widetilde{t}_{1}=0$
---
$\widetilde{\Delta t}_{1}$
---
$\widehat{\Delta t}_{1}$
---
$\widetilde{\Delta t}_{2}$
---
$\widehat{\Delta t}_{2}$
---
$\widetilde{\Delta t}_{3}$
---
$\widehat{\Delta t}_{3}$
---
$\widehat{\Delta t}_{N-1}$
---
$\widehat{\Delta t}_{N}$
---
$\widetilde{\Delta t}_{N}$
---
$\widetilde{\Delta t}_{N+1}$
---
$\widehat{t}_{1}$
---
$\widetilde{t}_{2}$
---
$\widehat{t}_{2}$
---
$\widetilde{t}_{3}$
---
$\widehat{t}_{3}$
---
$\widetilde{t}_{4}$
---
$\widehat{t}_{N-1}$
---
$\widetilde{t}_{N}$
---
$\widehat{t}_{N}$
---
$\widetilde{t}_{N+1}$
---
$\tau$
---
Figure 4: The entire analysis time $\tau$ is split into a number $N$ of
dissipation periods (indicated in blue) where a portion of the TPCL moves at
velocities of $\mathcal{O}\left(v_{\mathrm{cap}}\right)$, interspersed between
$N+1$ equilibrium stages where the entire TPCL moves at velocities of at most
$\mathcal{O}\left(v_{\mathrm{cv}}\right)$.

A number of mathematical constraints follow from the dynamic model of fluid
movement described above. Firstly, the analysis assumes that
$v_{\mathrm{cv}}\ll v_{\mathrm{cap}}$. By assuming that the fluids are
Newtonian with a stress defined by

$\mathsf{T}_{\mathrm{M}}=-p\mathsf{I}_{\mathrm{}}+\mu\left[\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}+(\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}})^{T}\right]$
(14)

an order of magnitude analysis on the augmented Navier-Stokes equations (as
defined by equations (1) and (2)) shows that capillary driven fluid velocities
can be limited by either viscous dissipation or momentum
acceleration/advection terms during any dissipation event. Hence, a
conservative estimate for the capillary driven velocity scale is given by
equating the interfacial stress term to either the viscous dissipation or
momentum acceleration/advection terms, giving

$v_{\mathrm{cap}}=\min\left(\frac{\sigma}{\mu},\sqrt{\frac{\sigma}{\rho
h_{\mathrm{rough}}}}\right)$ (15)

In this expression and subsequent order of magnitude analyses, properties such
as $\sigma$, $\mu$ and $\rho$ are order of magnitude estimates only, which for
most expressions can be taken as the maximum of the different phase properties
existing within the CV. Equation (15) places a constraint on the maximum
$v_{\mathrm{cv}}$ that can be used given that $v_{\mathrm{cv}}\ll
v_{\mathrm{cap}}$.

Interestingly, equation (15) predicts that capillary velocities are limited by
momentum acceleration/advection terms rather than viscous dissipation on most
practical surfaces. To illustrate, for a water droplet within air advancing
over a rough solid surface, using $\mu$ and $\rho$ from the water phase and
$\sigma$ as the surface tension coefficient between air and water gives
$\sigma/\mu\approx 72\,\mathrm{m/s}$. However, for all surface roughness
values $h_{\mathrm{rough}}\gtrapprox
h_{\mathrm{rough,crit}}=\mu^{2}/(\rho\sigma)=14\,\mathrm{nm}$ the capillary
velocity will be limited by the momentum acceleration/advection term
($\sqrt{\sigma/(\rho h_{\mathrm{rough}})}$) and hence will determine
$v_{\mathrm{cap}}$ for this system. Indeed, for a more typical surface
roughness of $h_{\mathrm{rough}}=10\,\mathrm{\mu m}$ equation (15) gives
$v_{\mathrm{cap}}\approx 3\,\mathrm{m/s}$ meaning that for practical surfaces
as long as $v_{\mathrm{cv}}$ is of the order cm/s or less $v_{\mathrm{cv}}\ll
v_{\mathrm{cap}}$ will be satisfied. We consider the variation of
$v_{\mathrm{cap}}$ in more detail in section 4.2. Note that even though the
local velocities existing during a dissipation event may be determined by a
balance between capillary and momentum acceleration/advection terms, the
Reynolds number for the motion near the surface roughness is not large ($\sim
20$ for the above $h_{\mathrm{rough}}=10\,\mathrm{\mu m}$ system). This means
that the size of the region where velocities are
$\mathcal{O}\left(v_{\mathrm{cap}}\right)$ during dissipation events is only
of $\mathcal{O}\left(h_{\mathrm{rough}}\right)$, and importantly these high
velocities will not exist at the circumference of the CV, located at roughly
$r_{\mathrm{cv}}$ from the TPCL.

A second constraint required by the dynamic model of fluid movement outlined
above is that $\tau_{\mathrm{cap}}\ll\tau$. To understand what conditions this
places on the physical system we first recognise that the time taken for each
dissipation event can be estimated from the interface velocity and distance
travelled during each event — $\widehat{\Delta
t}_{k}=\mathcal{O}\left(h_{\mathrm{rough}}/v_{\mathrm{cap}}\right)$ — and that
as there are $N$ dissipation events occurring during $\tau$,
$\tau_{\mathrm{cap}}=\mathcal{O}\left(Nh_{\mathrm{rough}}/v_{\mathrm{cap}}\right)$.
Recognising that the number of dissipation events occurring during the advance
duration is $\mathcal{O}\left(A_{\mathrm{cv}}/h_{\mathrm{rough}}^{2}\right)$
where $A_{\mathrm{cv}}=X_{\mathrm{cv}}l_{\mathrm{cv}}$, and that
$\tau=X_{\mathrm{cv}}/v_{\mathrm{cv}}$ leads to

$\frac{\tau_{\mathrm{cap}}}{\tau}=\mathcal{O}\left(\frac{l_{\mathrm{cv}}}{h_{\mathrm{rough}}}\frac{v_{\mathrm{cv}}}{v_{\mathrm{cap}}}\right)\ll
1$ (16)

Equation (16) can always be satisfied provided that $v_{\mathrm{cv}}$ is small
enough, which is obviously the case when determining $\theta_{\mathrm{a}}$ as
this angle is derived in the limit of $v_{\mathrm{cv}}\rightarrow 0$. If the
energy balance is being applied to moving interfaces however (i.e.,
$v_{\mathrm{cv}}\neq 0$) equation (16) places a constraint on the maximum
applicable $v_{\mathrm{cv}}$. This is discussed further in section 4.2.

Two final conceptual models concerning the continuous fluid movement occurring
during the equilibrium stages of the flow have to be developed in order to
apply the contact angle mechanical energy balance: These models are for
velocity gradient and pressure, both of which relate to the fluid stress.

For the velocity gradient, we note that during equilibrium stages the TPCL
advances at a speed of $\mathcal{O}\left(v_{\mathrm{cv}}\right)$ over the
solid, resulting in a velocity discontinuity at the solid surface if the
conventional continuum non-slip fluid boundary condition is applied. Indeed,
velocity profiles which satisfy the Navier-Stokes equations and that are
consistent with both a moving TPCL and the non-slip boundary condition are
available (Moffatt, 1964; Huh & Scriven, 1971), however these result in
velocity gradients near that TPCL that increase as $1/r$ (where $r$ is the
distance to the TPCL). We find that integrating these gradients over the
region surrounding the TPCL in our energy balance results in a energy
dissipation term for non-zero $v_{\mathrm{cv}}$ that diverges logarithmically
in an unphysical fashion, as others have found (Huh & Scriven, 1971).
Solutions proposed to this problem, which we invoke here, all involve removing
or limiting the stress (equivalently velocity gradients) within the fluid at
small distances ($h_{\mathrm{mol}}$) from the TPCL. Various justifications for
this limiting have been proposed (Huh & Scriven, 1971; Joanny & de Gennes,
1984; Petrov & Petrov, 1992), but most revolve around a breakdown of the
Newtonian or continuum model of a fluid at the TPCL where individual molecules
or particles within the fluid must ‘jump’ along the solid. The implications of
this limit are discussed further in section 4.2. Putting these concepts
together we hence estimate the velocity gradients existing within the fluid
during the equilibrium stages as

$\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}=\mathcal{O}\left(\frac{v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}\right)$
(17)

where as discussed $h_{\mathrm{mol}}$ is a small lengthscale related to the
molecular (or non-continuum) nature of the fluid.

For the continuous flow pressure variation, we perform an order of magnitude
on the single-phase Navier-Stokes equations, recognising that pressure
gradients may develop within either fluid in response to both viscous stress
and momentum acceleration/advection terms 222Gravitational effects are not
included when estimating the pressure variation as in general the fluid may
not be at rest (and hence may not experience a hydrostatic pressure
variation), however we note that inclusion of a gravitational contribution to
the pressure variation would not alter the final results of the energy balance
presented in section 3.3.10 as equivalent gravitational terms contribute to
this balance via the gravitational potential energy transport term
$\widetilde{T}_{3,k}$.. Given these assumptions we describe the pressure
variation within either fluid phase as

$\displaystyle\widetilde{\boldsymbol{\nabla}_{\mathrm{}}p}$
$\displaystyle=\mathcal{O}\left(\rho\boldsymbol{\nabla}_{\mathrm{}}\cdot\widetilde{\boldsymbol{v}_{\mathrm{}}}\widetilde{\boldsymbol{v}_{\mathrm{}}}\right)+\mathcal{O}\left(\mu\boldsymbol{\nabla}_{\mathrm{}}^{2}\widetilde{\boldsymbol{v}_{\mathrm{}}}\right)$
$\displaystyle=\mathcal{O}\left(\frac{\rho
v_{\mathrm{cv}}^{2}}{\max(r,h_{\mathrm{mol}})}\right)+\mathcal{O}\left(\frac{\mu
v_{\mathrm{cv}}}{\left[\max(r,h_{\mathrm{mol}})\right]^{2}}\right)$ (18)

where again the fluid stress has been limited within a distance of
$h_{\mathrm{mol}}$ from the TPCL. Noting that the equilibrium stage fluid
pressure $\widetilde{p}$ is relative to some point in the surrounding fluid
far from the TPCL, and that under these slow flow conditions there is a
potential pressure jump over the fluid interface due to the surrounding
interface curvature of $\mathcal{O}\left(\sigma/h_{\mathrm{surround}}\right)$,
we model the non-gravitational pressure variation within the CV during the
equilibrium stages as

$\widetilde{p}=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}\right)+\mathcal{O}\left(\frac{\mu
v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}\right)+\mathcal{O}\left(\frac{\sigma}{h_{\mathrm{surround}}}\right)+p_{0}$
(19)

where $p_{0}$ is some reference pressure located at a point away from the
TPCL, but within the vicinity of the CV. All models describing the conceptual
behaviour of the fluid phases during the equilibrium stages have now been
defined.

For the solid phase, to apply the contact angle mechanical energy balance we
need models that describe how the solid velocities and stresses vary within
this phase as the TPCL advances over the rough surface. In this study we
invoke the simplest possible model by assuming that velocities within the
solid are everywhere zero. With this assumption the energy balance becomes
independent of solid phase stresses. On physical grounds zero velocities can
be justified within the solid by assuming it is a yield-stress (or plastic)
material that does not experience a stress exceeding its yield-stress during
the analysis time. In reality this is probably justifiable for most solids
used in engineering applications, but for soft solids used in (eg) biomedical
applications energy dissipation within the solid phase may be significant.
This is certainly an area for future work that could be incorporated into the
presented mechanical energy balance framework but is not advanced here.

###  3.3 Applying the mechanical energy balance to the moving control volume

With the physical system defined and semi-quantitative models for how the
materials within the CV behave as it advances over the rough solid, we can now
apply the mechanical energy balance to find the advancing contact angle.

####  3.3.1 Formulating the contact angle mechanical energy balance

We start with deriving the most general form of the contact angle energy
balance by applying equation (10) to the moving CV and integrating it over a
time period from $t_{1}$ to $t_{2}$, giving

$T_{0}(t_{1},t_{2})=\sum_{i=1}^{6}T_{i}(t_{1},t_{2})$ (20)

where

$\displaystyle T_{0}(t_{1},t_{2})$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\left[E(t=t_{2})-E(t=t_{1})\right]$
(21) $\displaystyle E(t)$
$\displaystyle=\int_{V_{\mathrm{cv}}}\left(\frac{1}{2}\rho
v^{2}+\rho\hat{\Phi}\right)\,dV+\sum_{i<j}\sigma_{ij}A_{ij}$ (22)
$\displaystyle T_{1}(t_{1},t_{2})$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{S_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S},ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}\,dSdt$
(23) $\displaystyle T_{2}(t_{1},t_{2})$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\frac{1}{2}\rho
v^{2}(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}})\,dSdt$ (24)
$\displaystyle T_{3}(t_{1},t_{2})$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\rho\hat{\Phi}(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}})\,dSdt$
(25) $\displaystyle T_{4}(t_{1},t_{2})$
$\displaystyle=-\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{S_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S,}ij}\boldsymbol{n}_{\text{S,}ij}\boldsymbol{n}_{\text{S,}ij}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}\,dSdt$
(26) $\displaystyle T_{5}(t_{1},t_{2})$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{S_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}\,dSdt$
(27) $\displaystyle T_{6}(t_{1},t_{2})$
$\displaystyle=-\frac{1}{A_{\mathrm{cv}}}\int^{t_{2}}_{t_{1}}\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}\,dVdt$
(28)

In the derivation use has been made of
$\mathsf{I}_{\mathrm{}}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}=\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{}}$.

Next, defining the following notations that correspond to the entire, the
$k$th equilibrium and the $k$th dissipation stages of the advance as

$\displaystyle\overline{T}_{i}$ $\displaystyle=T_{i}(t_{1}=0,t_{2}=\tau)$ (29)
$\displaystyle\widetilde{T}_{i,k}$
$\displaystyle=T_{i}(t_{1}=\widetilde{t}_{k},t_{2}=\widetilde{t}_{k}+\widetilde{\Delta
t}_{k})$ (30) $\displaystyle\widehat{T}_{i,k}$
$\displaystyle=T_{i}(t_{1}=\widehat{t}_{k},t_{2}=\widehat{t}_{k}+{\widehat{\Delta
t}}_{k})$ (31)

respectively for $i=0$ to $6$, equation (20) is applied over the entire
advance period from time $0$ to $\tau$ to give

$\overline{T}_{0}=\sum_{i=1}^{6}\overline{T}_{i}$ (32)

Recognising that $\overline{T}_{i}$ for $i=1$ to $6$ are all integrals over
the total time period, using equation (13) these terms can be written as sums
of the corresponding terms from the equilibrium and dissipation stages, giving

$\overline{T}_{0}=\sum_{i=1}^{6}\left(\sum_{k=1}^{N+1}\widetilde{T}_{i,k}+\sum_{k=1}^{N}\widehat{T}_{i,k}\right)=\sum_{i=1}^{6}\sum_{k=1}^{N+1}\widetilde{T}_{i,k}+\sum_{k=1}^{N}\sum_{i=1}^{6}\widehat{T}_{i,k}$
(33)

To simplify the final term we apply the energy balance equation (20) to the
$k$th dissipation period, giving

$\widehat{T}_{0,k}=\sum_{i=1}^{6}\widehat{T}_{i,k}$ (34)

which when substituted back into equation (33) leads to

$\overline{T}_{0}=\sum_{i=1}^{6}\sum_{k=1}^{N+1}\widetilde{T}_{i,k}+\sum_{k=1}^{N}\widehat{T}_{0,k}$
(35)

This is the form of the energy balance that is used to evaluate the advancing
contact angle. It expresses the total change in mechanical energy within the
moving CV between the start and end of the advance as the sum of energy
transfers happening during each of the equilibrium stages plus the change in
mechanical energy existing within the CV that occurs over each of the
dissipation events.

We now examine each of the terms in equation (35) finding either their order
of magnitude, or for terms that prove to be significant, expressions that
allow their evaluation in terms of system properties.

####  3.3.2 Examining term $\overline{T}_{0}$

This term represents the change in mechanical energy within the CV between the
start and end of the advance. Defining the notation $\overline{\Delta
a}=a(t=\tau)-a(t=0)$, $\overline{T}_{0}$ can be expressed using equations (21)
and (22) as

$\displaystyle\overline{T}_{0}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\overline{\Delta E}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int_{V_{\mathrm{cv}}}\frac{1}{2}\overline{\Delta(\rho
v^{2})}dV+\frac{1}{A_{\mathrm{cv}}}\int_{V_{\mathrm{cv}}}\overline{\Delta\rho}\hat{\Phi}dV+\sum_{i<j}\frac{\sigma_{ij}\overline{\Delta
A_{ij}}}{A_{\mathrm{cv}}}$ (36)

The first term on the RHS (right hand side) of this equation captures changes
to the kinetic energy within the CV between the start and end of the advance.
As both $t=0$ and $t=\tau$ are within equilibrium periods, changes to $\rho
v^{2}$ within the CV between $t=0$ and $\tau$ will be limited to a fluid
volume that is within $\mathcal{O}\left(h_{\mathrm{rough}}\right)$ of the
solid surface. Further, as fluid velocities are
$\mathcal{O}\left(v_{\mathrm{cv}}\right)$ at both times, the first term can be
evaluated as $\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}r_{\mathrm{cv}}/X_{\mathrm{cv}}\right)$
after cancelling $l_{\mathrm{cv}}$. Similarly, for the second gravitational
potential energy term, changes to $\rho$ are also limited to the same volume
within $\mathcal{O}\left(h_{\mathrm{rough}}\right)$ of the solid. Recognising
that the gravitational potential function satisfies
$\hat{\Phi}=\boldsymbol{g}_{\mathrm{}}\cdot\boldsymbol{x}_{\mathrm{}}$, a
maximum magnitude for this potential function within the CV is
$\hat{\Phi}=\mathcal{O}\left(gl_{\mathrm{cv}}\right)$ where
$g=|\boldsymbol{g}_{\mathrm{}}|$ is the gravitational constant. Finally, for
the third term in equation (36), which represents the change in surface
potential energy within the CV between the start and end of the advance, the
change in area of each interface will be of
$\mathcal{O}\left(h_{\mathrm{rough}}l_{\mathrm{cv}}\right)$, being composed of
changes to interfacial areas that occur around the TPCL, as well as changes to
the average of each $A_{ij}$ associated with solid interfaces under each of
the fluid phases due to the (possibly) random nature of the surface defects.

With these assumptions the change in mechanical energy over the duration of
the advance is evaluated as

$\displaystyle\overline{T}_{0}$ $\displaystyle=\mathcal{O}\left(\frac{\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}h_{\mathrm{rough}}}{X_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\rho
gr_{\mathrm{cv}}l_{\mathrm{cv}}h_{\mathrm{rough}}}{X_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\sigma
h_{\mathrm{rough}}}{X_{\mathrm{cv}}}\right)$
$\displaystyle=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
gl_{\mathrm{cv}}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\frac{\sigma
h_{\mathrm{rough}}}{X_{\mathrm{cv}}}\right)$ (37)

where in the last line we have used
$\mathcal{O}\left(r_{\mathrm{cv}}\right)=\mathcal{O}\left(X_{\mathrm{cv}}\right)$.
The above three terms are later compared to other terms present in equation
(35) to gauge their significance in determining the advancing contact angle
$\theta_{\mathrm{a}}$.

####  3.3.3 Examining term $\sum_{k=1}^{N}\widehat{T}_{0,k}$

In a similar fashion to $\overline{T}_{0}$, but here applied over each
dissipation period, $\widehat{T}_{0,k}$ represents the change in mechanical
energy within the CV occurring over the period of the $k$th dissipation event.
Defining $\widehat{\Delta a}_{k}=a(t=\widehat{t}_{k}+{\widehat{\Delta
t}}_{k})-a(t=\widehat{t}_{k})$ this change in energy can be expressed as

$\displaystyle\widehat{T}_{0,k}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\widehat{\Delta E}_{k}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\int_{V_{\mathrm{cv}}}\frac{1}{2}\widehat{\Delta(\rho
v^{2})}_{k}dV+\frac{1}{A_{\mathrm{cv}}}\int_{V_{\mathrm{cv}}}\widehat{\Delta\rho}_{k}\hat{\Phi}dV+\sum_{i<j}\frac{\sigma_{ij}\widehat{\Delta
A_{ij}}_{k}}{A_{\mathrm{cv}}}$ (38)

where $\widehat{\Delta A_{ij}}_{k}$ is the change in area of interface $ij$
that is contained within the CV and that occurs over the $k$th dissipation
event.

In order to evaluate the magnitude of the terms appearing in equation (38), we
return to our conceptual model for how the fluid and interfaces behave during
dissipation events. During a dissipation event, an area of the TPCL ‘depins’
from a particular surface defect and moves at a capillary-driven velocity to a
new ‘equilibrium’ interface position. These dissipation events cause
interfacial areas to change by
$\mathcal{O}\left(h_{\mathrm{rough}}^{2}\right)$, and as interfacial
curvatures resulting from the roughness extend by
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$ into the fluid, cause fluid
properties to change within a volume of
$\mathcal{O}\left(h_{\mathrm{rough}}^{3}\right)$ near the TPCL. Hence changes
to the kinetic energy and density within a volume of
$\mathcal{O}\left(h_{\mathrm{rough}}^{3}\right)$ caused by each dissipation
event will contribute to the first two terms on the RHS of equation (38),
while changes of $\mathcal{O}\left(h_{\mathrm{rough}}^{2}\right)$ to the
interfacial areas of each phase combination due to each dissipation event will
contribute to the third term on the RHS of this equation.

Concurrently, over the duration of each dissipation event (${\widehat{\Delta
t}}_{k}$) continuous movement of the TPCL still occurs across the rough
surface, and this movement also contributes to the terms on the RHS of
equation (38). Specifically, within
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$ of the TPCL, or a volume of
$\mathcal{O}\left({\widehat{\Delta
t}}_{k}v_{\mathrm{cv}}l_{\mathrm{cv}}h_{\mathrm{rough}}\right)$, there will be
a change in kinetic energy and density of the fluid occurring due to the
continuous TPCL movement which will add contributions to the first two terms
on the RHS of equation (38). Similarly, for the third term on the RHS of
equation (38), there will also be a change in interfacial areas of
$\mathcal{O}\left({\widehat{\Delta
t}}_{k}v_{\mathrm{cv}}l_{\mathrm{cv}}\right)$ due to the continuous TPCL
movement that also needs to be included.

Hence, summing changes due to both the specific dissipation event and
continuous TPCL movement occurring during each dissipation period, the
magnitude of $\widehat{T}_{0,k}$ can be expressed as

$\displaystyle\widehat{T}_{0,k}$ $\displaystyle=\mathcal{O}\left(\frac{\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}^{3}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\rho
v_{\mathrm{cv}}^{3}h_{\mathrm{rough}}{\widehat{\Delta
t}}_{k}}{X_{\mathrm{cv}}}\right)$ $\displaystyle+\mathcal{O}\left(\frac{\rho
gh_{\mathrm{rough}}^{4}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\rho
gl_{\mathrm{cv}}v_{\mathrm{cv}}h_{\mathrm{rough}}{\widehat{\Delta
t}}_{k}}{X_{\mathrm{cv}}}\right)$ $\displaystyle+\mathcal{O}\left(\frac{\sigma
h_{\mathrm{rough}}^{2}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\sigma
v_{\mathrm{cv}}{\widehat{\Delta t}}_{k}}{X_{\mathrm{cv}}}\right)$ (39)

where each pair of the six terms corresponds to the first, second and third
terms appearing on the RHS of equation (38), respectively, with the first of
each pair corresponding to changes caused by the specific capillary-driven
dissipation event, and the second of each pair corresponding to the changes
due to continuous TPCL movement that occurs during each dissipation event
period.

Returning to equation (35), it is actually the sum of $\widehat{T}_{0,k}$ from
all $N$ dissipation events that is required in the contact angle mechanical
energy balance. Performing this sum on equation (39), while noting
$\sum_{k=1}^{N}{\widehat{\Delta t}}_{k}=\tau_{\mathrm{cap}}$,
$N=\mathcal{O}\left(X_{\mathrm{cv}}l_{\mathrm{cv}}/h_{\mathrm{rough}}^{2}\right)$
and $v_{\mathrm{cv}}\tau=X_{\mathrm{cv}}$ leads to

$\displaystyle\sum_{k=1}^{N}\widehat{T}_{0,k}$
$\displaystyle=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\frac{\tau_{\mathrm{cap}}}{\tau}\right)$
$\displaystyle+\mathcal{O}\left(\rho
gh_{\mathrm{rough}}^{2}\right)+\mathcal{O}\left(\rho
gl_{\mathrm{cv}}h_{\mathrm{rough}}\frac{\tau_{\mathrm{cap}}}{\tau}\right)$
$\displaystyle+\mathcal{O}\left(\sigma\right)+\mathcal{O}\left(\sigma\frac{\tau_{\mathrm{cap}}}{\tau}\right)$
(40)

Noting from equation (16) that
$\mathcal{O}\left(\tau_{\mathrm{cap}}/\tau\right)\ll 1$, terms involving this
ratio can be neglected in comparison to other terms, and recognising that the
$\mathcal{O}\left(\sigma\right)$ term in the above originated from the final
term of equation (38), we arrive at

$\sum_{k=1}^{N}\widehat{T}_{0,k}=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
gh_{\mathrm{rough}}^{2}\right)+\widehat{\Delta\sigma}$ (41)

Here the specific dissipation event surface energy change per area traversed
has been defined as

$\widehat{\Delta\sigma}=\sum_{k=1}^{N}\sum_{i<j}\sigma_{ij}\frac{\widehat{\Delta
A_{ij}}_{k}}{A_{\mathrm{cv}}}$ (42)

This variable represents the sum of changes to potential surface energies
occurring within the CV due to all capillary-driven dissipation events, per
projected area of solid traversed. This variable has a magnitude of
$\mathcal{O}\left(\sigma\right)$. It is similar to the $W_{d}$ term used in
Joanny & de Gennes (1984) and will become key in determining
$\theta_{\mathrm{a}}$ from the contact angle energy analysis.

The next six terms all correspond to energy transfers that occur during the
equilibrium stages of the advance.

####  3.3.4 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{1,k}$

This term represents transport of surface potential energy through the
boundary of the CV during the equilibrium stages due only to movement of the
CV. Introducing the shorthand notation $\int_{\widetilde{\Delta
t}_{k}}=\int_{\widetilde{t}_{k}}^{\widetilde{t}_{k}+\widetilde{\Delta
t}_{k}}$, applying equation (23) over $\widetilde{\Delta t}_{k}$ gives

$\sum_{k=1}^{N+1}\widetilde{T}_{1,k}=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S},ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt$
(43)

Further, noting that
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}$ on the ends of
the CV ($S_{\mathrm{cv,end}}$), and that $\sigma_{ij}$ is non-zero only within
three thin regions on the circumference of the CV ($S_{\mathrm{cv,cir}}$)
where the interfaces cross its boundary, equation (43) can be written as the
sum of three terms

$\displaystyle\sum_{k=1}^{N+1}\widetilde{T}_{1,k}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,top}}}\sigma_{12}\delta_{\text{S},12}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt$
$\displaystyle+\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,bl}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S},ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt$
$\displaystyle+\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,br}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S},ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt$
(44)

Here, as indicated in Figure 5, $S_{\mathrm{cv,top}}$, $S_{\mathrm{cv,bl}}$
and $S_{\mathrm{cv,br}}$ are three thin regions located on
$S_{\mathrm{cv,cir}}$ that just contain the fluid interface (top), interfaces
associated with the solid surface under phase $\bigcirc$1 (bl = bottom left)
and interfaces associated with the solid surface under phase $\bigcirc$2 (br =
bottom right), respectively.

$S_{\mathrm{cv,br}}$
---
$S_{\mathrm{cv,bl}}$
---
$S_{\mathrm{cv,top}}$
---
$S_{\mathrm{cv,end}}$
---
$S_{\mathrm{cv,cir}}$
---
$S_{\mathrm{cv,end}}$
---
$\theta_{\mathrm{a}}$
---
$\bigcirc$1
---
$\bigcirc$2
---
$\bigcirc$S
---
$\boldsymbol{v}_{\mathrm{cv}}$
---
Figure 5: Three thin regions on the circumference of the CV
($S_{\mathrm{cv,cir}}$) are defined that contain all of the phase interfaces
that intersect with this boundary.

For the first term on the RHS of equation (44), as
$\boldsymbol{n}_{\mathrm{cv}}=-\cos\theta_{\mathrm{a}}\boldsymbol{e}_{1}+\sin\theta_{\mathrm{a}}\boldsymbol{e}_{3}$
where the fluid interface crosses the CV boundary at $S_{\mathrm{cv,top}}$,
and as $\boldsymbol{v}_{\mathrm{cv}}=v_{\mathrm{cv}}\boldsymbol{e}_{1}$,
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}=-v_{\mathrm{cv}}\cos\theta_{\mathrm{a}}$
at this location. Further, following the interface description outlined in
Section 3.1,
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{n}_{\mathrm{S},ij}=0$ at this
location at all times, and as $S_{\mathrm{cv,top}}$ is infinitely thin in the
circumferential direction at this location equation (4) shows that
$\int_{S_{\mathrm{cv,top}}}\delta_{\text{S,}ij}dS=l_{\mathrm{cv}}$. Hence the
first term on the RHS of equation (44) becomes

$\displaystyle\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,top}}}\sigma_{12}\delta_{\text{S,}ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt$
$\displaystyle=-\frac{\sigma_{12}v_{\mathrm{cv}}\cos\theta_{\mathrm{a}}}{X_{\mathrm{cv}}}\sum_{k=1}^{N+1}\widetilde{\Delta
t}_{k}$
$\displaystyle=-\left(1-\frac{\tau_{\mathrm{cap}}}{\tau}\right)\sigma_{12}\cos\theta_{\mathrm{a}}$
(45)

The second and third terms on the RHS of equation (44) represent the transport
of surface potential energy associated with the rough solid surface out of,
and into, respectively, the control volume as it advances across the solid.
Note that these surface energies may include contributions from micro-droplets
or bubbles confined within the surface roughness, as well as from the surface
energy between the solid and adjacent bulk fluid phase.

Focusing our explanation on the second term on the RHS of (44) that involves
an integral over $S_{\mathrm{cv,bl}}$, we first note that as
$S_{\mathrm{cv,bl}}$ only has to include phase interfaces that are associated
with the surface roughness, the circumferential height of $S_{\mathrm{cv,bl}}$
is only of $\mathcal{O}\left(h_{\mathrm{rough}}\right)$. Further as
$h_{\mathrm{rough}}\ll r_{\mathrm{cv}}$ (via equation (11)) we have
$\boldsymbol{n}_{\mathrm{cv}}=-\boldsymbol{e}_{1}$ at this location and
consequently
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}=-v_{\mathrm{cv}}$.
Splitting the temporal integration using equation (13) leads to

$\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,bl}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S,}ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt=\\\
-\frac{v_{\mathrm{cv}}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\int_{0}^{\tau}\int_{S_{\mathrm{cv,bl}}}\delta_{\text{S,}ij}dSdt\\\
-\frac{v_{\mathrm{cv}}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\sum_{k=1}^{N+1}\int_{{\widehat{\Delta
t}}_{k}}\int_{S_{\mathrm{cv,bl}}}\delta_{\text{S},ij}dSdt$ (46)

For the first term on the RHS of this equation, we define a new volume
$\overleftarrow{V}$ that is created by sweeping $S_{\mathrm{cv,bl}}$ over the
solid surface for the duration of the advance. Formally we define a coordinate
$x_{\mathrm{cv}}=tv_{\mathrm{cv}}$ that increases in the direction of
$\boldsymbol{e}_{1}$ and measures the progress of the CV as it moves over the
solid surface, giving

$\displaystyle
v_{\mathrm{cv}}\sum_{i<j}\sigma_{ij}\int_{0}^{\tau}\int_{S_{\mathrm{cv,bl}}}\delta_{\text{S,}ij}dSdt$
$\displaystyle=\sum_{i<j}\sigma_{ij}\int_{0}^{X_{\mathrm{cv}}}\int_{S_{\mathrm{cv,bl}}}\delta_{\text{S,}ij}dSdx_{\mathrm{cv}}$
$\displaystyle=\sum_{i<j}\sigma_{ij}\int_{\overleftarrow{V}}\delta_{\text{S,}ij}dV$
$\displaystyle=\sum_{i<j}\sigma_{ij}\overleftarrow{A}_{ij}$ (47)

where $\overleftarrow{A}_{ij}$ is the area of each interface type $ij$
associated with the rough solid surface that leaves the CV during its advance
over the solid.

For the second term on the RHS of equation (46) we calculate its order of
magnitude rather than derive an expression, noting that
$\sum_{k=1}^{N}{\widehat{\Delta t}}_{k}=\tau_{\mathrm{cap}}$ and
$\int_{S_{\mathrm{cv,bl}}}\delta_{\text{S},ij}dS=\mathcal{O}\left(l_{\mathrm{cv}}\right)$.
Incorporating these expressions and equation (47) back into equation (46)
gives

$\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,bl}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S,}ij}\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{v}_{\mathrm{cv}}dSdt=-\sum_{i<j}\sigma_{ij}\frac{\overleftarrow{A}_{ij}}{A_{\mathrm{cv}}}+\mathcal{O}\left(\frac{\tau_{\mathrm{cap}}}{\tau}\sigma\right)$
(48)

completing the evaluation of the second term on the RHS of equation (44).

Finally, performing an analogous calculation for the third term on the RHS of
equation (44) for the solid surface interfaces $\overrightarrow{A}_{ij}$ that
lie under phase $\bigcirc$2 and that enter the CV through $S_{\mathrm{cv,br}}$
during the advance, and noting again as per equation (16) that terms involving
the factor $\tau_{\mathrm{cap}}/\tau$ can be neglected relative to those that
do not include this term, the surface potential energy transport term from
equation (44) can be written as

$\sum_{k=1}^{N+1}\widetilde{T}_{1,k}=-\sigma_{12}\cos\theta_{\mathrm{a}}+\overrightarrow{\sigma}-\overleftarrow{\sigma}$
(49)

where the compound surface energies associated with the solid interface that
is leaving (behind and to the left of the TPCL) and entering (in front and to
the right of the TPCL) the advancing CV are defined by

$\overleftarrow{\sigma}=\sum_{i<j}\sigma_{ij}\frac{\overleftarrow{A}_{ij}}{A_{\mathrm{cv}}}\quad\text{and}\quad\overrightarrow{\sigma}=\sum_{i<j}\sigma_{ij}\frac{\overrightarrow{A}_{ij}}{A_{\mathrm{cv}}},$
(50)

respectively. We note that these compound surface energies incorporate any
micro-bubbles or droplets that could possibly be contained within the solid
surface roughness. For the leaving energy ($\overleftarrow{\sigma}$) these
could be formed during the discussed dissipation events that occur as the TPCL
sweeps over rough solid. For the entering energy ($\overrightarrow{\sigma}$)
these could be pre-existing within the surface roughness, possibly a result of
previous wetting processes (such as receding back over a previously wetted
surface).

####  3.3.5 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{2,k}$

This term represents the transport of kinetic energy into the CV during the
equilibrium stages, and is given by

$\sum_{k=1}^{N+1}\widetilde{T}_{2,k}=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\frac{1}{2}\rho
v^{2}(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}})dSdt$ (51)

Noting from section 3.2 that velocities are zero within the solid phase and
within the equilibrium stages within the fluid are
$\widetilde{\boldsymbol{v}_{\mathrm{}}}=\mathcal{O}\left(v_{\mathrm{cv}}\right)$,
and that the fluid areas of $S_{\mathrm{cv,cir}}$ and $S_{\mathrm{cv,end}}$
have areas of $\mathcal{O}\left(r_{\mathrm{cv}}l_{\mathrm{cv}}\right)$ and
$\mathcal{O}\left(r_{\mathrm{cv}}^{2}\right)$, respectively, the magnitude of
this term is given by

$\sum_{k=1}^{N+1}\widetilde{T}_{2,k}=\mathcal{O}\left((\tau-\tau_{\mathrm{cap}})\rho
v_{\mathrm{cv}}^{3}\frac{(r_{\mathrm{cv}}^{2}+l_{\mathrm{cv}}r_{\mathrm{cv}})}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\right)=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}\right)$ (52)

In deriving the final term in this expression we have used
$\mathcal{O}\left(l_{\mathrm{cv}}\right)=\mathcal{O}\left(r_{\mathrm{cv}}\right)$
and neglected a term containing $\tau_{\mathrm{cap}}/\tau$ relative to one
that does not, consistent with equation (16).

####  3.3.6 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{3,k}$

This term represents the transport of gravitational potential energy into the
CV during the equilibrium stages, and is given by

$\sum_{k=1}^{N+1}\widetilde{T}_{3,k}=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv}}}\boldsymbol{n}_{\mathrm{cv}}\cdot\rho\hat{\Phi}(\boldsymbol{v}_{\mathrm{cv}}-\boldsymbol{v}_{\mathrm{}})dSdt$
(53)

In a very similar fashion to the last term this term has a magnitude given by

$\sum_{k=1}^{N+1}\widetilde{T}_{3,k}=\mathcal{O}\left((\tau-\tau_{\mathrm{cap}})\rho
v_{\mathrm{cv}}g\frac{r_{\mathrm{cv}}^{2}}{X_{\mathrm{cv}}}\right)=\mathcal{O}\left(\rho
gr_{\mathrm{cv}}^{2}\right)$ (54)

where we have additionally assumed that
$\hat{\Phi}=\mathcal{O}\left(gr_{\mathrm{cv}}\right)$ over
$S_{\mathrm{cv,cir}}$ and
$\hat{\Phi}=\mathcal{O}\left(gl_{\mathrm{cv}}\right)$ over
$S_{\mathrm{cv,end}}$.

####  3.3.7 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{4,k}$

The next term is a contributor to the work that the interfaces outside the CV
do on the material inside the CV. Using equations (26) and (30) this term is
given by

$\sum_{k=1}^{N+1}\widetilde{T}_{4,k}=-\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv}}}\sum_{i<j}\sigma_{ij}\delta_{\text{S,}ij}\boldsymbol{n}_{\text{S,}ij}\boldsymbol{n}_{\text{S,}ij}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dSdt$
(55)

As the integral contains the surface delta function $\delta_{\text{S,}ij}$,
like $\widetilde{T}_{1,k}$, only four component surfaces of $S_{\mathrm{cv}}$
give non-zero contributions to the integral: Namely $S_{\mathrm{cv,top}}$,
$S_{\mathrm{cv,bl}}$, $S_{\mathrm{cv,br}}$ and $S_{\mathrm{cv,end}}$. We
evaluate
$\boldsymbol{n}_{\text{S,}ij}\boldsymbol{n}_{\text{S,}ij}:\widetilde{\boldsymbol{v}_{\mathrm{}}}\boldsymbol{n}_{\mathrm{cv}}$
on each of these surfaces.

On $S_{\mathrm{cv,top}}$,
$\boldsymbol{n}_{\mathrm{cv}}\cdot\boldsymbol{n}_{\text{S,}ij}=0$ as here the
fluid interface is perfectly flat and normal to the CV boundary (as discussed
in section 3.1) so this surface makes no contribution to equation (55). For
$S_{\mathrm{cv,bl}}$ and $S_{\mathrm{cv,br}}$ the interfaces associated with
these surface areas have $\boldsymbol{v}_{\mathrm{}}=0$ (the solid velocity)
at the CV circumference, so these surfaces also make no contribution to
$\widetilde{T}_{4,k}$. Finally, over $S_{\mathrm{cv,end}}$, at distances from
the TPCL that are much larger than
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$,
$\boldsymbol{n}_{\text{S,}ij}\cdot\boldsymbol{n}_{\mathrm{cv}}=0$ as like over
$S_{\mathrm{cv,top}}$ here the fluid interface is flat and normal to the CV
boundary. However at distances from the TPCL on $S_{\mathrm{cv,end}}$ that are
of $\mathcal{O}\left(h_{\mathrm{rough}}\right)$,
$\boldsymbol{n}_{\text{S,}ij}\cdot\boldsymbol{n}_{\mathrm{cv}}$ is not zero as
in this region the interface has curvature of
$\mathcal{O}\left(1/h_{\mathrm{rough}}\right)$, as described in section 3.1.
Further, within this region the fluid velocity
$\widetilde{\boldsymbol{v}_{\mathrm{}}}=\mathcal{O}\left(v_{\mathrm{cv}}\right)$
as the TPCL may be advancing over the solid here during these equilibrium
stages. Hence there is a contribution to the integral in equation (55) from
the inner part of $S_{\mathrm{cv,end}}$ only, giving overall

$\sum_{k=1}^{N+1}\widetilde{T}_{4,k}=\mathcal{O}\left((\tau-\tau_{\mathrm{cap}})\sigma
v_{\mathrm{cv}}\frac{h_{\mathrm{rough}}}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\right)=\mathcal{O}\left(\sigma\frac{h_{\mathrm{rough}}}{l_{\mathrm{cv}}}\right)$
(56)

where as per equation (16) a term involving the factor
$\tau_{\mathrm{cap}}/\tau$ has been neglected.

####  3.3.8 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{5,k}$

This term represents the work that the material stresses
$\mathsf{T}_{\mathrm{M}}$ are doing on the material inside the CV. Using
equations (27) and (30) this term is defined as

$\sum_{k=1}^{N+1}\widetilde{T}_{5,k}=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dSdt$
(57)

Noting that within the solid
$\boldsymbol{v}_{\mathrm{}}=\boldsymbol{0}_{\mathrm{}}$ (as discussed section
3.2), and applying the fluid Newtonian stress equation (14) this term can be
written as

$\displaystyle\sum_{k=1}^{N+1}\widetilde{T}_{5,k}$
$\displaystyle=-\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,fluid}}}\widetilde{p}\widetilde{\boldsymbol{v}_{\mathrm{}}}\cdot\boldsymbol{n}_{\mathrm{cv}}dSdt$
$\displaystyle+\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{S_{\mathrm{cv,fluid}}}\mu\left[\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}+(\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}})^{T}\right]:\widetilde{\boldsymbol{v}_{\mathrm{}}}\boldsymbol{n}_{\mathrm{cv}}dSdt$
(58)

where we have used the identity
$\mathsf{I}_{\mathrm{}}:\widetilde{\boldsymbol{v}_{\mathrm{}}}\boldsymbol{n}_{\mathrm{cv}}=\widetilde{\boldsymbol{v}_{\mathrm{}}}\cdot\boldsymbol{n}_{\mathrm{cv}}$
and $S_{\mathrm{cv,fluid}}$ represents the surface of the CV within the fluid
phase. Equilibrium stage fluid properties $\widetilde{p}$ and
$\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}$ are
relevant during these time intervals and are substituted from the steady-state
order of magnitude expressions of equations (19) and (17), respectively, as
discussed in section 3.2, giving

$\sum_{k=1}^{N+1}\widetilde{T}_{5,k}=-\frac{(\tau-\tau_{\mathrm{cap}})}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\times\\\
\int_{S_{\mathrm{cv,fluid}}}\left[\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}+\frac{\mu
v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}+\frac{\sigma}{h_{\mathrm{surround}}}\right)+p_{0}\right]\widetilde{\boldsymbol{v}_{\mathrm{}}}\cdot\boldsymbol{n}_{\mathrm{cv}}dS\\\
+\frac{(\tau-\tau_{\mathrm{cap}})}{X_{\mathrm{cv}}l_{\mathrm{cv}}}\int_{S_{\mathrm{cv,fluid}}}\mathcal{O}\left(\frac{\mu
v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}\right):\widetilde{\boldsymbol{v}_{\mathrm{}}}\boldsymbol{n}_{\mathrm{cv}}dS$
(59)

For the term involving the reference pressure $p_{0}$, as this is constant at
any given time it can come out of the integral, leaving
$\widetilde{\boldsymbol{v}_{\mathrm{}}}\cdot\boldsymbol{n}_{\mathrm{cv}}$
which is zero when integrated over $S_{\mathrm{cv,fluid}}$ as the fluid is
incompressible (using Gauss’ theorem). Hence this reference pressure term does
not contribute to $\widetilde{T}_{5,k}$. For the remainder we gather like
terms and evaluate in an order of magnitude sense over the circumference and
ends of the CV that are within the fluid region as

$\displaystyle\sum_{k=1}^{N+1}\widetilde{T}_{5,k}$
$\displaystyle=\mathcal{O}\left\\{(1-\frac{\tau_{\mathrm{cap}}}{\tau})\frac{1}{l_{\mathrm{cv}}}\left[\int_{S_{\mathrm{cv,cir}}}\left(\rho
v_{\mathrm{cv}}^{2}+\frac{\mu
v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}+\frac{\sigma}{h_{\mathrm{surround}}}\right)dS\right.\right.$
$\displaystyle\quad\quad\left.\left.+\int_{S_{\mathrm{cv,end}}}\left(\rho
v_{\mathrm{cv}}^{2}+\frac{\mu
v_{\mathrm{cv}}}{\max(r,h_{\mathrm{mol}})}+\frac{\sigma}{h_{\mathrm{surround}}}\right)dS\right]\right\\}$
$\displaystyle=\mathcal{O}\left\\{(1-\frac{\tau_{\mathrm{cap}}}{\tau})\frac{1}{l_{\mathrm{cv}}}\left[\left(\rho
v_{\mathrm{cv}}^{2}+\frac{\mu
v_{\mathrm{cv}}}{r_{\mathrm{cv}}}+\frac{\sigma}{h_{\mathrm{surround}}}\right)r_{\mathrm{cv}}l_{\mathrm{cv}}\right.\right.$
$\displaystyle\quad\quad\left.\left.+\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}^{2}+\int_{0}^{h_{\mathrm{mol}}}\frac{\mu
v_{\mathrm{cv}}}{h_{\mathrm{mol}}}rdr+\int_{h_{\mathrm{mol}}}^{r_{\mathrm{cv}}}\frac{\mu
v_{\mathrm{cv}}}{r}rdr+\frac{\sigma
r_{\mathrm{cv}}^{2}}{h_{\mathrm{surround}}}\right]\right\\}$
$\displaystyle=\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}\right)+\mathcal{O}\left(\mu
v_{\mathrm{cv}}\right)+\mathcal{O}\left(\frac{\sigma
r_{\mathrm{cv}}}{h_{\mathrm{surround}}}\right)$ (60)

where for the last line we have used
$\mathcal{O}\left(r_{\mathrm{cv}}\right)=\mathcal{O}\left(l_{\mathrm{cv}}\right)$
and employed equations (11) and (16) to neglect comparatively small terms.

####  3.3.9 Examining term $\sum_{k=1}^{N+1}\widetilde{T}_{6,k}$

This final term required for the evaluation of equation (35) represents the
rate of dissipation occurring within the CV during the equilibrium periods,
and is evaluated in a very similar manner to $\widetilde{T}_{5,k}$. Using
equations (28) and (30) this term is defined as

$\sum_{k=1}^{N+1}\widetilde{T}_{6,k}=-\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{V_{\mathrm{cv}}}\mathsf{T}_{\mathrm{M}}:\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}\boldsymbol{n}_{\mathrm{cv}}dVdt$
(61)

Noting again that within the solid
$\boldsymbol{v}_{\mathrm{}}=\boldsymbol{0}_{\mathrm{}}$, and applying the
fluid Newtonian stress equation (14) we find

$\displaystyle\sum_{k=1}^{N+1}\widetilde{T}_{6,k}$
$\displaystyle=\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{V_{\mathrm{cv,fluid}}}\widetilde{p}\mathsf{I}_{\mathrm{}}:\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}dVdt$
$\displaystyle-\frac{1}{A_{\mathrm{cv}}}\sum_{k=1}^{N+1}\int_{\widetilde{\Delta
t}_{k}}\int_{V_{\mathrm{cv,fluid}}}\mu\left[\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}+(\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}})^{T}\right]:\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}dVdt$
(62)

For the pressure term we use
$\mathsf{I}_{\mathrm{}}:\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}=\boldsymbol{\nabla}_{\mathrm{}}\cdot\widetilde{\boldsymbol{v}_{\mathrm{}}}=0$
as the fluid is incompressible, removing this entire integral. Substituting
the equilibrium stage fluid velocity gradient from equation (17) into the
remaining dissipation integral and evaluating in an order of magnitude sense
leads to

$\displaystyle\sum_{k=1}^{N+1}\widetilde{T}_{6,k}$
$\displaystyle=\mathcal{O}\left(\frac{\mu(\tau-\tau_{\mathrm{cap}})}{A_{\mathrm{cv}}}\int_{V_{\mathrm{cv,fluid}}}\frac{v_{\mathrm{cv}}^{2}}{\left[\max(r,h_{\mathrm{mol}})\right]^{2}}dV\right)$
$\displaystyle=\mathcal{O}\left[\mu
v_{\mathrm{cv}}\left(1-\frac{\tau_{\mathrm{cap}}}{\tau}\right)\left(\int_{0}^{h_{\mathrm{mol}}}\frac{r}{h_{\mathrm{mol}}^{2}}dr+\int_{h_{\mathrm{mol}}}^{r_{\mathrm{cv}}}\frac{1}{r}dr\right)\right]$
$\displaystyle=\mathcal{O}\left[\mu
v_{\mathrm{cv}}\ln\left(\frac{r_{\mathrm{cv}}}{h_{\mathrm{mol}}}\right)\right]$
(63)

where again relatively small terms have been neglected via equations (11) and
(16).

####  3.3.10 Final contact angle mechanical energy balance

With all relevant terms defined, we substitute equations (37), (41), (49),
(52), (54), (56), (60), (63) into equation (35) giving

$\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
gl_{\mathrm{cv}}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\frac{\sigma
h_{\mathrm{rough}}}{X_{\mathrm{cv}}}\right)=\\\ \mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
gh_{\mathrm{rough}}^{2}\right)+\widehat{\Delta\sigma}-\sigma_{12}\cos\theta_{\mathrm{a}}+\overrightarrow{\sigma}-\overleftarrow{\sigma}\\\
+\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}\right)+\mathcal{O}\left(\rho
gr_{\mathrm{cv}}^{2}\right)+\mathcal{O}\left(\sigma\frac{h_{\mathrm{rough}}}{l_{\mathrm{cv}}}\right)\\\
+\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}\right)+\mathcal{O}\left(\mu
v_{\mathrm{cv}}\right)+\mathcal{O}\left(\frac{\sigma
r_{\mathrm{cv}}}{h_{\mathrm{surround}}}\right)+\mathcal{O}\left[\mu
v_{\mathrm{cv}}\ln\left(\frac{r_{\mathrm{cv}}}{h_{\mathrm{mol}}}\right)\right]$
(64)

Gathering like terms and neglecting any terms that are relatively small due to
the separation of lengthscales equation (11) yields

$\widehat{\Delta\sigma}-\sigma_{12}\cos\theta_{\mathrm{a}}+\overrightarrow{\sigma}-\overleftarrow{\sigma}=\\\
+\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}\right)+\mathcal{O}\left(\rho
gr_{\mathrm{cv}}^{2}\right)+\mathcal{O}\left[\mu
v_{\mathrm{cv}}\ln\left(\frac{r_{\mathrm{cv}}}{h_{\mathrm{mol}}}\right)\right]$
(65)

Equivalently this contact angle energy balance can be expressed as

$\sigma_{12}\cos\theta_{\mathrm{a}}=\widehat{\Delta\sigma}+\overrightarrow{\sigma}-\overleftarrow{\sigma}$
(66)

under conditions where the following inequalities hold

$\mathcal{O}\left(\frac{\rho
v_{\mathrm{cv}}^{2}r_{\mathrm{cv}}}{\sigma}\right),\mathcal{O}\left(\frac{\rho
gr_{\mathrm{cv}}^{2}}{\sigma}\right),\mathcal{O}\left[\frac{\mu
v_{\mathrm{cv}}}{\sigma}\ln\left(\frac{r_{\mathrm{cv}}}{h_{\mathrm{mol}}}\right)\right]\ll
1$ (67)

####  3.3.11 Interpretation of specific dissipation event surface energy
change $\widehat{\Delta\sigma}$

The contact angle energy balance given in equation (66) involves the specific
dissipation event surface energy change $\widehat{\Delta\sigma}$, which as
previously defined is the sum of all surface energy changes occurring within
the CV as a result of individual dissipation events (equation (42)).
Consistent with the physical model of how the contact line advances over the
rough solid as outlined in section 3.2, surface energy is dissipated to heat
(via viscous stresses) during each dissipation event resulting in a negative
$\widehat{\Delta\sigma}$. Hence an alternative nomenclature for this term that
is consistent with this contact line movement model (as discussed in the
introduction, (e.g. Joanny & de Gennes, 1984)) is to define the specific
dissipation per projected area of solid traversed due to surface roughness as
$D\approx-\widehat{\Delta\sigma}$. In this section we formalise this
relationship by performing an order of magnitude analysis on an energy balance
conducted over all dissipation events occurring during the advance period to
show that the surface energy liberated via $\widehat{\Delta\sigma}$ is in fact
dissipated to heat during these periods.

In order to perform this dissipation event analysis, we need models for how
the materials within the CV behave during the dissipation periods. As for the
equilibrium periods, we assume that the solid does not deform (following the
justifications from section 3.2) giving $\boldsymbol{v}_{\mathrm{}}=0$ within
the solid during these times. For the fluids we require a description of the
velocity $\widehat{\boldsymbol{v}_{\mathrm{}}}$, velocity gradient
$\widehat{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}$ and
pressure $\widehat{p}$ existing during these times. For the velocity, noting
that the continuous equilibrium velocities still exist within the dissipation
periods, and assuming that significant capillary driven velocities only exist
within an region of size $\mathcal{O}\left(h_{\mathrm{rough}}\right)$ local to
each dissipation event (as the Reynolds number of these motions is not large),
we define

$\widehat{\boldsymbol{v}_{\mathrm{}}}(\boldsymbol{x}_{\mathrm{}})=\begin{cases}\widetilde{\boldsymbol{v}_{\mathrm{}}}+\mathcal{O}\left(v_{\mathrm{cap}}\right)&\text{if
}r_{\mathrm{cap}}<\mathcal{O}\left(h_{\mathrm{rough}}\right)\\\
\widetilde{\boldsymbol{v}_{\mathrm{}}}&\text{otherwise}\end{cases}$ (68)

where
$r_{\mathrm{cap}}=|\boldsymbol{x}_{\mathrm{}}-\boldsymbol{x}_{\mathrm{cap}}|$
is the distance to the centre of the relevant $k$th dissipation event, and
$\boldsymbol{x}_{\mathrm{cap}}$ is the location of the particular $k$th
dissipation event that is centred on the centreline of the CV 333For
notational simplicity here and in subsequent dissipation event model
definitions we do not indicate what specific dissipation event variables such
as $\widehat{\boldsymbol{v}_{\mathrm{}}}$, $\boldsymbol{x}_{\mathrm{cap}}$ and
$r_{\mathrm{cap}}$ refer to. For the velocity gradient we similarly define

$\widehat{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}(\boldsymbol{x}_{\mathrm{}})=\begin{cases}\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}+\mathcal{O}\left(\frac{v_{\mathrm{cap}}}{\max(r_{\mathrm{cap}},h_{\mathrm{mol}})}\right)&\text{if
}r_{\mathrm{cap}}<\mathcal{O}\left(h_{\mathrm{rough}}\right)\\\
\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}&\text{otherwise}\end{cases}$
(69)

where as per the equilibrium velocity gradient model of equation (17) we limit
the stress generated at the moving TPCL that is within
$\mathcal{O}\left(h_{\mathrm{mol}}\right)$ of the solid. Finally for pressure
we perform another order of magnitude analysis on the Navier-Stokes equations,
but now recognise that capillary induced pressure changes occur over the
interface that is deforming over length
$\mathcal{O}\left(h_{\mathrm{rough}}\right)$ during these times, giving

$\widehat{p}(\boldsymbol{x}_{\mathrm{}})=\begin{cases}\widetilde{p}+\mathcal{O}\left(\rho
v_{\mathrm{cap}}^{2}\right)+\mathcal{O}\left(\frac{\mu
v_{\mathrm{cap}}}{\max(r_{\mathrm{cap}},h_{\mathrm{mol}})}\right)+\mathcal{O}\left(\frac{\sigma}{h_{\mathrm{rough}}}\right)&\text{if
}r_{\mathrm{cap}}<\mathcal{O}\left(h_{\mathrm{rough}}\right)\\\
\widetilde{p}&\text{otherwise}\end{cases}$ (70)

This completes the dissipation event material specifications.

The dissipation event energy analysis now largely mirrors that conducted for
the equilibrium stages (as detailed in sections 3.3.4 to 3.3.9) by summing the
individual dissipation event energy balances of equation (34) over all $N$
dissipation events. The details of this analysis are contained within Appendix
B, with the final conclusion of the energy balance being

$\displaystyle\widehat{\Delta\sigma}$
$\displaystyle=\mathcal{O}\left(\frac{\tau_{\mathrm{cap}}}{\tau}\sum_{i=1}^{6}\sum_{k=1}^{N+1}\widetilde{T}_{i,k}\right)+\mathcal{O}\left(\rho
v_{\mathrm{cv}}^{2}h_{\mathrm{rough}}\right)+\mathcal{O}\left(\rho
gh_{\mathrm{rough}}^{2}\right)$ $\displaystyle+\mathcal{O}\left(\frac{\rho
v_{\mathrm{cap}}^{2}h_{\mathrm{rough}}^{2}}{l_{\mathrm{cv}}}\right)+\mathcal{O}\left(\frac{\rho
gh_{\mathrm{rough}}^{3}}{l_{\mathrm{cv}}}\right)+\mathcal{O}\left(\sigma\frac{h_{\mathrm{rough}}}{l_{\mathrm{cv}}}\right)$
$\displaystyle+\mathcal{O}\left(\mu
v_{\mathrm{cap}}\frac{h_{\mathrm{rough}}}{l_{\mathrm{cv}}}\right)+\mathcal{O}\left[\mu
v_{\mathrm{cap}}\ln\left(\frac{h_{\mathrm{rough}}}{h_{\mathrm{mol}}}\right)\right]$
(71)

Under conditions for which the inequalities of equation (67) hold, the
equilibrium terms represented by $\sum_{k=1}^{N+1}\widetilde{T}_{i,k}$ in the
above have a maximum magnitude of $\mathcal{O}\left(\sigma\right)$, which is
the same as that of $\widehat{\Delta\sigma}$. As
$\mathcal{O}\left(\frac{\tau_{\mathrm{cap}}}{\tau}\right)\ll 1$ via equation
(16), it follows that the equilibrium terms make no significant contribution
to equation (71) and can be neglected. For the remaining terms, using the
lengthscales assumption equation (11) in combination with equation (67) yields
after simplifications

$\widehat{\Delta\sigma}=\mathcal{O}\left(\frac{\rho
v_{\mathrm{cap}}^{2}h_{\mathrm{rough}}^{2}}{l_{\mathrm{cv}}}\right)+\mathcal{O}\left[\mu
v_{\mathrm{cap}}\ln\left(\frac{h_{\mathrm{rough}}}{h_{\mathrm{mol}}}\right)\right]$
(72)

Tracing back through the dissipation energy analysis presented in Appendix B
shows that the first and second terms on the RHS of equation (72) represent
the transport of kinetic energy (and associated pressure work) through the
ends of the CV due to dissipation event velocities, and the viscous
dissipation of energy that is converted to heat around the moving TPCL during
dissipation events, respectively. The ratio of these two terms depends on the
roughness scale capillary velocity Reynolds number as well as various
lengthscale ratios, so that in general neither can assumed to be dominant. The
second term directly represents energy conversion to heat during the capillary
driven dissipation events, as envisaged. While the first term does not
directly relate to energy dissipation, it does represent the transport of
energy through the ends of the CV and hence parallel to the TPCL. As the
analysis shows that no significant energy transport occurs over the
circumference of the CV during the dissipation events, the energy represented
by the first term in equation (72) remains within the vicinity of the contact
line, and so will eventually be dissipated to heat via the second term in the
same equation, but in an adjacent control volume. Hence, overall equation (72)
shows that the energy liberated by surface changes that occur during each
dissipation event is dissipated to heat within the vicinity of the TPCL, and
we are justified in using the nomenclature $D=-\widehat{\Delta\sigma}$ to
represent the specific energy dissipation per area travelled due to surface
roughness.

It should be reiterated however that the discussion in this section does not
change the way in which $D$ (or equivalently $-\widehat{\Delta\sigma}$) is
evaluated, being defined solely in terms of surface energy changes that occur
over all dissipation events. Indeed a strength of this framework is that $D$
can be calculated based solely on these surface energy changes, rather than
needing to know the dynamic and transient details of each individual
dissipation event.

##  4 Discussion

###  4.1 Summary of key results

The contact angle mechanical energy balance equations are summarised in table
1. The requirements of inequality equation (67) have been combined with
previous lengthscale assumptions and interpreted as conditions on
$r_{\mathrm{cv}}$ and $v_{\mathrm{cv}}$.
$\mathcal{O}\left(l_{\mathrm{cv}}\right)=\mathcal{O}\left(r_{\mathrm{cv}}\right)$
is applied throughout. The limiting CV size variable $r_{\mathrm{cv,grav}}$ is
the capillary length that in this analysis has originated from limiting the
transport of gravitational potential energy during the equilibrium stages of
the energy balance. The three limiting control volume velocities
$v_{\mathrm{cv,ke}}$, $v_{\mathrm{cv,vis}}$ and $v_{\mathrm{cv,cap}}$ have
originated from limiting the transport of kinetic energy through the CV
boundaries during the equilibrium stages, limiting the rate of viscous
dissipation occurring at a molecular level around the moving TPCL during the
equilibrium stages, and ensuring that the duration of dissipation events
occurring within the CV is small compared to the duration of interface
movement, respectively. As per section (3.3.11), the specific change in
surface potential energy occurring during the dissipation periods has been
defined in terms of $D$, being the specific dissipation due to roughness.

$\sigma_{12}\cos\theta_{\mathrm{a}}=\overrightarrow{\sigma}-\overleftarrow{\sigma}-D$
(73)

subject to

$\displaystyle h_{\mathrm{mol}}\ll h_{\mathrm{rough}}\ll r_{\mathrm{cv}}\ll
h_{\mathrm{surround}},r_{\mathrm{cv,grav}}$ (74) $\displaystyle
v_{\mathrm{cv}}\ll\min\left(v_{\mathrm{cv,ke}},v_{\mathrm{cv,vis}},v_{\mathrm{cv,cap}}\right)$
(75)

where

$\displaystyle\overleftarrow{\sigma}=\sum_{i<j}\sigma_{ij}\frac{\overleftarrow{A}_{ij}}{A_{\mathrm{cv}}}\quad\text{and}\quad\overrightarrow{\sigma}=\sum_{i<j}\sigma_{ij}\frac{\overrightarrow{A}_{ij}}{A_{\mathrm{cv}}}$
(50) $\displaystyle
D=-\sum_{k=1}^{N}\sum_{i<j}\sigma_{ij}\frac{\widehat{\Delta
A_{ij}}_{k}}{A_{\mathrm{cv}}}$ (76) $\displaystyle
v_{\mathrm{cap}}=\min\left(\frac{\sigma}{\mu},\sqrt{\frac{\sigma}{\rho
h_{\mathrm{rough}}}}\right)$ (15) $\displaystyle
v_{\mathrm{cv,ke}}=\sqrt{\frac{\sigma}{\rho r_{\mathrm{cv}}}}$ (77)
$\displaystyle
v_{\mathrm{cv,vis}}=\sqrt{\frac{\sigma}{\mu\ln\left(r_{\mathrm{cv}}/h_{\mathrm{mol}}\right)}}$
(78) $\displaystyle
v_{\mathrm{cv,cap}}=v_{\mathrm{cap}}\frac{h_{\mathrm{rough}}}{r_{\mathrm{cv}}}$
(79) $\displaystyle r_{\mathrm{cv,grav}}=\sqrt{\frac{\sigma}{\rho g}}$ (80)

Table 1: Equation summary for the contact angle mechanical energy balance.
$D$, $\overleftarrow{\sigma}$, $\overrightarrow{\sigma}$ are the specific
dissipation due to roughness, leaving compound solid surface energy and
entering compound solid surface energy, respectively.

It is illustrative to consider some simple examples of the application of
equation (73). In each of these cases phase $\bigcirc$1 is advancing to the
right at $\theta_{\mathrm{a}}$ while displacing phase $\bigcirc$2 (as per
Figure 2):

  * •

Youngs: For a flat surface the leaving compound surface energy is
$\overleftarrow{\sigma}=\sigma_{\mathrm{1s}}$, the entering compound surface
energy is $\overrightarrow{\sigma}=\sigma_{\mathrm{2s}}$ and the specific
roughness dissipation rate $D=0$ as no dissipation events occur during the
TPCL advance. Application of equation (73) gives
$\sigma_{12}\cos\theta_{\mathrm{a}}=\sigma_{\mathrm{2s}}-\sigma_{\mathrm{1s}}$,
or via Young’s equation, $\theta_{\mathrm{a}}=\theta_{\mathrm{e}}$.

  * •

Wenzel: Here the leaving surface energy is
$\overleftarrow{\sigma}=r\sigma_{\mathrm{1s}}$ and the entering surface energy
is $\overrightarrow{\sigma}=r\sigma_{\mathrm{2s}}$, where $r$ is the roughness
or specific solid area (total solid area per projected solid area).
Application of equation (73) gives
$\cos\theta_{\mathrm{a}}=r\cos\theta_{\mathrm{e}}-D/\sigma_{12}$ which is
consistent with Wenzel’s equation, but augmented by a roughness dissipation
term. Here increasing the roughness $r$ can either increase or decrease the
advancing angle depending on $\theta_{\mathrm{e}}$, however increasing the
roughness dissipation $D$ always increases the contact angle.

  * •

Cassie: In this general case $f_{1}$ and $f_{2}$ are defined as the wetted and
non-wetted liquid areas under the droplet per projected area, respectively,
with $r$ defined as per the Wenzel case. Under these conditions the leaving
surface energy is
$\overleftarrow{\sigma}=-f_{1}\sigma_{\mathrm{12}}\cos\theta_{\mathrm{e}}+f_{2}\sigma_{\mathrm{12}}+r\sigma_{\mathrm{2s}}$
and the entering surface energy is
$\overleftarrow{\sigma}=r\sigma_{\mathrm{2s}}$. Substituting into equation
(73) gives
$\cos\theta_{\mathrm{a}}=f_{1}\cos\theta_{\mathrm{e}}-f_{2}-D/\sigma_{12}$
which is consistent with Cassie’s analysis(Cassie & Baxter, 1944), but
augmented with a dissipation term. For a droplet resting on the top of flat
structures (for example ‘Fakir’ droplets existing on photolithography based
surfaces) this equation can be further simplified using $f_{2}=1-f_{1}$.

In general the specific roughness dissipation $D$ and leaving compound surface
energy $\overrightarrow{\sigma}$ cannot be guessed from the solid topology
alone and instead must be modelled from an understanding of the fluid dynamics
occurring within the CV, or measured. The entering compound surface energy
$\overleftarrow{\sigma}$ depends on the history of the surface and must
similarly be modelled or measured. Hence the mechanical energy balance
equation (73) is best viewed as a framework within which results from dynamic
interfacial modelling or measurements can be incorporated for specific
combinations of solid topologies and fluid phases. These activities are the
topic of on-going work in our lab.

In terms of calculating the receding angle, as per Figure 2 when phase
$\bigcirc$1 is advancing over the solid, phase $\bigcirc$2 is receding. Hence
substituting $\theta_{\mathrm{a}}=\pi-\theta_{\mathrm{r}}$ into equation (73)
and letting the compound surface energies under phases $\bigcirc$1 and
$\bigcirc$2 while phase $\bigcirc$2 is receding equal
$\overrightarrow{\sigma}_{\mathrm{r}}=\overleftarrow{\sigma}$ and
$\overleftarrow{\sigma}_{\mathrm{r}}=\overrightarrow{\sigma}$, respectively,
yields
$\sigma_{12}\cos\theta_{\mathrm{r}}=\overrightarrow{\sigma}_{\mathrm{r}}-\overleftarrow{\sigma}_{\mathrm{r}}+D_{\mathrm{r}}$,
with $D_{\mathrm{r}}$ now referring to the dissipation occurring while phase
$\bigcirc$2 recedes. As expected an increase in $D_{\mathrm{r}}$ leads to a
decrease in $\theta_{\mathrm{r}}$.

###  4.2 Range of validity

One strength of this energy balance approach for predicting contact angles is
that it not only provides a framework for predicting the angle, but also
specifies under what conditions the predicted angle will be valid.
Specifically, for equation (73) to be representative of a real system all
conditions specified by equations (74) and (75) must be satisfied. Continuing
the discussion from section 3.2, we examine what physical limitations these
conditions place on applying the theory to the common water/air system, with
results shown in Figure 6.

In terms of lengthscales, evaluating $r_{\mathrm{cv,grav}}$ from equation (80)
requires that $r_{\mathrm{cv}}\ll r_{\mathrm{cv,grav}}=2.7\,\mathrm{mm}$,
which is approximately satisfied if we adopt
$r_{\mathrm{cv}}=r_{\mathrm{cv,max}}=r_{\mathrm{cv,grav}}/10=0.27\,\mathrm{mm}$
for the water system (that is, one order of magnitude less). From the
separation of lengthscales equation (74) this places an upper constraint on
the roughness applicable under the theory of $h_{\mathrm{rough}}\lessapprox
27\,\mathrm{\mu m}$, which is indicated by the right bound of the shaded area
in Figure 6. Hence, gravitational effects place an upper limit on the size of
roughness applicable under this theory. At the small scales we also require
that $h_{\mathrm{mol}}\ll h_{\mathrm{rough}}$, and given that the molecular
size of water is around $0.27\,\mathrm{nm}$ this places a lower limit on
$h_{\mathrm{rough}}$ of around $2.7\,\mathrm{nm}$ as indicated by the left
bound of the shaded area in Figure 6. Note that some theories use a larger
length cut-off of $100\,\mathrm{nm}$ in their analysis citing the influence of
van der Waals and double-layer forces at these lengthscales(de Gennes, 1985):
While still continuum forces, these are not accounted for in our
hydrodynamically-based framework. Non-continuum effects such as thermal
fluctuations that could influence interface topologies at very small
lengthscales are also not considered. Hence, while the lower roughness size
limit of this theory requires further validation, overall the conclusion is
that CAH can be significant even for surfaces that have very small sized
roughness (at least in an engineering sense), a conclusion that is supported
by experiments(Delmas et al., 2011; Fetzer & Ralston, 2011).

It is interesting that aside from these lengthscale and velocity constraints
(discussed below), there is no reference to the absolute size of the surface
roughness contained within the energy analysis. While models of $D$,
$\overleftarrow{\sigma}$ and $\overrightarrow{\sigma}$ may in some cases
depend on the absolute size of the surface roughness, it is likely that for
many fluid combinations and solid topologies the advancing and receding
contact angles do not depend on the absolute size of the roughness, but only
on the topology of the roughness. As discussed in the introduction, this is
consistent with a growing number of observations (Öner & McCarthy, 2000;
Dorrer & Rühe, 2008; Li et al., 2016; Jiang et al., 2019).

$v_{\mathrm{cv}}$(m/s)

$h_{\mathrm{rough}}$(m)$v_{\mathrm{cap}}$$v_{\mathrm{cv,cap}}$$v_{\mathrm{cv,ke}}$$v_{\mathrm{cv,vis}}$$v_{\mathrm{cv,cap^{\prime}}}$$v_{\mathrm{cv,ke^{\prime}}}$$v_{\mathrm{cv,vis^{\prime}}}$$r_{\mathrm{cv,grav}}$$r_{\mathrm{cv,max}}$$h_{\mathrm{mol}}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{-10}$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$
Figure 6: Regions of validity (shaded) for the contact angle energy analysis
based on a typical water/air system. The velocity limits
$v_{\mathrm{cv,cap}}$, $v_{\mathrm{cv,ke}}$ and $v_{\mathrm{cv,vis}}$ and the
entire shaded validity region (dark grey plus light grey) assume that the CV
size depends on the solid roughness according to
$r_{\mathrm{cv}}/h_{\mathrm{rough}}=10$. Conversely, the velocity limits
$v_{\mathrm{cv,cap^{\prime}}}$, $v_{\mathrm{cv,ke^{\prime}}}$ and
$v_{\mathrm{cv,vis^{\prime}}}$ and light grey validity region assume that the
CV size is constant at the maximum allowable value of
$r_{\mathrm{cv}}=r_{\mathrm{cv,max}}=r_{\mathrm{cv,grav}}/10$.

In terms of the speed that the CV can travel at for the theory still to be
valid (that is, $v_{\mathrm{cv}}$), the constraints partly depend on how the
size of the CV is defined. We consider two scenarios: In the first the CV size
decreases with $h_{\mathrm{rough}}$ such that
$r_{\mathrm{cv}}/h_{\mathrm{rough}}=10$; In the second the CV size is kept
constant at the maximum value of
$r_{\mathrm{cv}}=r_{\mathrm{cv,max}}=0.27\,\mathrm{mm}$ as determined by
equation (74). The various limiting velocities corresponding to the first
decreasing CV size case are indicated by $v_{\mathrm{cv,cap}}$,
$v_{\mathrm{cv,ke}}$ and $v_{\mathrm{cv,vis}}$ in Figure 6, while the limiting
velocities corresponding to the second constant CV size case are indicated by
$v_{\mathrm{cv,cap^{\prime}}}$, $v_{\mathrm{cv,ke^{\prime}}}$ and
$v_{\mathrm{cv,vis^{\prime}}}$. Similarly, the ranges of theory validity are
shown in the figure by the entire shaded region for the first case, and by the
light shaded region for the second case. In practice the CV size determines at
what distance from the solid (or TPCL) we are measuring the ‘macroscopic’
advancing angle, $\theta_{\mathrm{a}}$.

For both CV size cases considered the capillary velocity constraint, which is
an interpretation of equation (16) and requires that the total dissipation
event time during the advance ($\tau_{\mathrm{cap}}$) is less than the total
advance time ($\tau$), is the only constraint that limits the advance velocity
over the entire range of applicable roughness scales. As was discussed in
section 3.2 and shown in Figure 6, the capillary velocity $v_{\mathrm{cap}}$
for this water system increases as $h_{\mathrm{rough}}$ decreases, reaching
$\sim 72\,\mathrm{m/s}$ below
$h_{\mathrm{rough}}=h_{\mathrm{rough,crit}}\approx 14\,\mathrm{nm}$. For the
first decreasing CV size case this velocity constraint allows
$v_{\mathrm{cv}}$ to remain at practically large values (of $\sim
0.72\,\mathrm{m/s}$) even at the lowest limit of applicable
$h_{\mathrm{rough}}$. For the second constant CV size case however the
capillary velocity constraint decreases the viable $v_{\mathrm{cv}}$ as
$h_{\mathrm{rough}}$ decreases, reaching $\sim 70\,\mathrm{\mu m/s}$ at the
lowest applicable $h_{\mathrm{rough}}\approx 2.7\,\mathrm{nm}$.

It is notable that for most roughness sizes and under either CV size scenarios
the velocity of the interface under which the theory is valid, and hence
predicts that the advancing angle will be equal to the static advancing angle
of the system, is in practical terms large. At the highest valid solid
roughnesses that are in the tens of micron range the advancing angle will
remain as the static angle until the interface is moving faster than several
cm/s. At the lowest valid roughnesses the maximum valid advance velocity
depends on the choice of CV size, but even for the more restrictive larger CV
size case the maximum advance velocity is still around $0.1\,\mathrm{mm/s}$ at
the smallest roughnesses, which is certainly an experimentally accessible
value for performing contact angle measurements. In general these results
emphasise that contact angle hysteresis depends more on surface roughness than
interface velocity, for most practical surfaces experiencing interface
movement for most engineering applications.

It is also interesting that the maximum allowable interface speed is in
general a function of the CV size, or put differently, that the predicted
dynamic contact angle is a function of the distance from the TPCL (and hence
solid) that the angle is measured. This has been noted previously in the
context of dynamic contact angle models that are based on the concept of
viscous dissipation (de Gennes, 1985). As previously discussed, for the second
CV size regime as specified by the light grey area in Figure 6, the upper
velocity constraint is determined by $v_{\mathrm{cv,cap^{\prime}}}$, which can
alternatively be interpreted as a constraint on the maximum number of solid
defects existing per solid area over which the CV advances
($A_{\mathrm{cv}}$). This constraint is represented by equation (16). However,
it is possible that this constraint could be relaxed if a different approach
to the energy balance was adopted: Rather than requiring that the total time
for dissipation events be relatively small, an alternative approach that was
advanced in earlier versions of this theory (Harvie et al., 2008) is to locate
dissipation events not only in time but also in space within the CV, allowing
more than one dissipation event to be occurring within the analysis
concurrently. This approach is mathematically more onerous to produce rigorous
conditions on theory validity, but is worthy of future research as it could
relax the $v_{\mathrm{cv,cap}}$ constraint on $v_{\mathrm{cv}}$ under
conditions where the CV size is of a macroscopic (i.e. mm) size but the solid
roughness much smaller. Note that if the $v_{\mathrm{cv,cap}}$ restraint were
relaxed, the next constraint on the dynamic contact angle would come from
$v_{\mathrm{cv,ke}}$ which is related to kinetic energy transport through the
CV volume, rather than $v_{\mathrm{cv,vis}}$ which is concerned with viscous
dissipation occurring within the CV volume (and more specifically around the
TPCL). Interestingly previous theories have focused on predicting dynamic
contact angles from knowledge of viscous or molecular energy dissipation
occurring around the TPCL(Moffatt, 1964; Huh & Scriven, 1971; Voinov, 1976; de
Gennes, 1985; Petrov & Petrov, 1992), rather than considering kinetic energy
transport within the TPCL region — in contrast to this framework that suggests
that kinetic energy transport is more important than viscous dissipation in
predicting dynamic contact angles. Validating this result and producing new
models of dynamic contact angle is another area for future research suggested
by this energy conservation framework.

###  4.3 Other model limitations

In addition to the advance velocity, roughness size and CV size constraints
discussed in the previous section, other assumptions used in the derivation of
the energy balance theory will not be valid for some systems, but could form
avenues for future extensions of the theory.

Inherent in the contact angle framework, and enabling multiphase mechanical
energy balance presented in section 2, is the assumption that there is no
irreversible work involved in interface creation (or destruction). Physically
irreversibilities could be present in a system via chemical or molecular
effects associated with interface formation, including effects due to
surfactants or other types of surface active molecules. In the context of
dynamic contact angle modelling, some work has been done in incorporating
molecular dissipation and adsorption when predicting advancing contact
angles(Blake & Haynes, 1969; de Gennes, 1985; Brochard-Wyart & Gennes, 1992;
Mohammad Karim, 2022). Following this work and recognising the parallels
between these interfacial formation irreversibilities and dissipation due to
TPCL movement, it is likely that this irreversible interfacial work will
appear in the contact angle energy balance as a dissipation term that is in
addition to the hydodynamically based $D$, but this does require reformulation
of the general mechanical energy balance and application to the moving contact
line CV to show this rigorously. A further assumption used by the model is
that the solid is chemically homogeneous. This assumption could be relaxed in
a straight-forward manner by including more than one solid interface type in
the contact angle energy analysis.

Another limitation of the presented work is that fluids are assumed to be
incompressible, and no dissolution or evaporation of the fluids is permitted.
Both of these assumptions are physically limiting when one of the phases is a
gas. Physically a gas will behave differently to a liquid in cases where a
micro-bubble is formed within the solid roughness, caused by the specific
dynamics of the fluid interface as it advances over solid defects. Within a
formed micro-bubble the laplace pressure will be high
($\mathcal{O}\left(\sigma/h_{\mathrm{rough}}\right)$) which will cause the
density of gas to increase and hence the volume of the micro-bubble to
decrease. How significant these affects are depends on the bubble size,
interfacial tension and gas equation of state. Further, high pressures within
the bubble will drive gas dissolution into the surrounding fluid phase,
further reducing the volume of the formed micro-bubbles. These changes to the
micro-bubble size will in turn affect the leaving compound solid surface
energy $\overleftarrow{\sigma}$, affecting the predicted advancing contact
angle. In terms of the energy framework, fluid compressibility could be
implemented quite easily by relaxing the
$\boldsymbol{\nabla}_{\mathrm{}}\cdot\boldsymbol{v}_{\mathrm{}}=0$ constraint
used when evaluating the boundary pressure work terms. Gas dissolution is more
complex however, as this process is transient, and so the framework would need
to recognise over what timescale the dissolution process is taking place, and
hence where in the energy balances the changes in (particularly) interfacial
energy should occur. Micro-droplet evaporation is a related phenomena not
accounted for by the theory that could also change $\overleftarrow{\sigma}$
and $\overrightarrow{\sigma}$ depending on the wetting history of the solid
surface and volatility of the fluids. It would be interesting to account for
dissolution and evaporation effects within this framework as both introduce a
transient or history effect to the wetting behaviour of the system. Again,
this suggests avenues for future theoretical and experimental research.

A further limitation of the theory relates to the size of roughness
considered. As presented the roughness is characterised by a single
lengthscale $h_{\mathrm{rough}}$ which represents both the height of surface
defects and the spacing between them. An extension to this theory could
distinguish between these two lengthscales, allowing better analysis of
surfaces composed of a dilute number of strong defects, or surfaces that are
gently undulating. A second assumption relating to the size of the defects is
that only one lengthscale is considered in the analysis, whereas real surfaces
(particularly biological inspired super-hydrophobic surfaces) can possess a
hierarchical range of roughnesses(Feng et al., 2002). A potential extension of
the theory would be to use a cascade of CV sizes to predict the contact angles
at each lengthscale, with the equilibrium contact angles used at each
lengthscale determined from contact angle predictions performed using smaller
lengthscale CVs. In this way the macroscopic contact angle could be determined
from knowledge of roughness topologies at each lengthscale.

Finally the stress model for the solid material used in this study effectively
implies that there are no energy changes occurring within the CV that are
related to the solid phase. More complex solid stress models could certainly
be incorporated into the framework, representing (for example) the advance of
fluids over soft or semi-liquid materials. Introducing deformable solids would
not only introduce additional dissipation, but would also require additional
dynamical modelling or measuring of the fluid, solid and interfacial
behaviours occurring around the TPCL. This is also the topic of on-going work.

##  5 Conclusion

Starting from a statement of momentum conservation, a mechanical energy
conservation framework has been derived that allows the contact angle
hysteresis (CAH) range to be predicted from knowledge of the interfacial
dynamics that occur around an advancing three phase contact line (TPCL).
Unlike most previous works the analysis is not specific to a particular
wetting regime (e.g., Cassie or Wenzel) or particular surface structure (e.g.,
holes, poles, periodic, dilute etc). The analysis also resolves a number of
questions about wetting on rough surfaces that have been the source of
confusion in the literature:

  * •

A suitable control volume (CV) can be defined that is anchored to the TPCL
around which mechanical energy conservation can be rigorously performed. This
shows that at least for the static CAH range, and conditional upon the CV
lengthscale and speed satisfying four conditions (i.e., equations (77-80)),
the advancing and receding angles are independent of the surrounding flow.
This is significant as it means that the energy of the entire flow system
(e.g., the entire drop) need not be considered when predicting CAH, resolving
extensive discussion within the field(Gao & McCarthy, 2007; McHale, 2007; Gao
& McCarthy, 2007; Nosonovsky, 2007; Panchagnula & Vedantam, 2007; Marmur,
2022).

  * •

By adopting semi-quantitative models for how the fluids and phase interfaces
behave within the advancing CV that are based on experimental observation, we
have shown from energy conservation that CAH directly results from the
‘pinning’ and ‘depinning’ or ‘jumping’ of the fluid interface, deriving a
general expression for this dissipation in terms of surface areas pre and post
each jumping event. No dissipation or CAH results from contact line distortion
(in isolation), showing (in the language of Joanny & de Gennes (1984)) that
only ‘strong’ defects (as opposed to ‘weak’ defects) cause CAH. This is
consistent with the initial observations of Joanny & de Gennes (1984) but
conflicts with some subsequent studies (Pomeau & Vannimenus, 1985; Robbins &
Joanny, 1987; Öpik, 2000).

  * •

Studies that consider how a fluid interface advances over a rough surface have
previously inferred what the macroscopic CAH range is from the results, but
without clear justification(Cox, 1983). The present work shows how interfacial
and energy analyses can be related: Namely, interfacial dynamic models can be
used to predict dissipation and generated (leaving) compound solid surface
energies which can then be fed back into the derived energy conservation
framework to predict the CAH range. Only for specific simple systems will the
CAH range inferred from interfacial analysis match that predicted via energy
conservation and observed at the macroscopic scale.

  * •

While not the central focus of this study, the energy conservation analysis
suggests that when calculating dynamic contact angles kinetic energy transport
around the TPCL is at least as significant as hydrodynamic energy dissipation
occurring around the TPCL. This requires more theoretical and experimental
confirmation, but could explain some of the discrepancies between existing
dynamic contact angle models and experimental observation (Mohammad Karim,
2022).

Many limitations of the framework have also been discussed and could form the
basis of future extensions, including considering compressible, evaporating or
dissolving fluids, irreversible work associated with interface creation or
destruction, solid chemical heterogeneity, soft solids, and non-negligible CV
advance speeds. These are topics of on-going work in our group.

##  6 Nomenclature

Arabic:
---
$A_{ij}$ |  area of $ij$ interface type
$\overrightarrow{A}_{ij}$ |  area of $ij$ interface type entering the control volume from the right
$\overleftarrow{A}_{ij}$ |  area of $ij$ interface type leaving the control volume on the left
$A_{\mathrm{cv}}$ |  projected area of solid surface swept by CV over advance duration
$\widehat{\Delta A_{ij}}_{k}$ |  change in interfacial area $ij$ during dissipation period $k$
$D$ |  total energy dissipation occurring during advance
$\boldsymbol{e}_{k}$ |  coordinate unit vector in direction $k$
$\boldsymbol{g}_{\mathrm{}}$ |  gravity vector
$h_{\mathrm{rough}}$ |  lengthscale of solid surface roughness
$h_{\mathrm{rough,crit}}$ |  lengthscale of solid surface roughness that delimits viscous and inertial fluid flow regimes
$h_{\mathrm{surround}}$ |  lengthscale of surround flow
$h_{\mathrm{mol}}$ |  lengthscale of molecular or non-continuum effects within the fluid
$\mathsf{I}_{\mathrm{}}$ |  identity tensor
$l_{\mathrm{cv}}$ |  length of CV
$\boldsymbol{n}_{\mathrm{cv}}$ |  outward normal to surface of CV
$\boldsymbol{n}_{\text{S,}ij}$ |  unit normal vector of $ij$ interface type, directed into phase $i$
$m$ |  number of material phases present within control volume
$p$ |  pressure
$p_{0}$ |  reference pressure
$\widehat{p}$ |  pressure within the fluid during a dissipation period
$\widetilde{p}$ |  pressure within the fluid during an equilibrium period
$r_{\mathrm{cv}}$ |  radius of CV
$r_{\mathrm{cap}}$ |  distance to centre of dissipation event
$r_{\mathrm{cv,grav}}$ |  capillary length
$r_{\mathrm{cv,max}}$ |  maximum applicable CV size
$S_{\mathrm{cv}}$ |  surface of CV
$S_{\mathrm{cv,cir}}$ |  circumferential surface of CV
$S_{\mathrm{cv,end}}$ |  end surface of CV
$S_{\mathrm{cv,br}}$ |  circumferential surface of CV at the bottom right containing compound solid interfaces
$S_{\mathrm{cv,bl}}$ |  circumferential surface of CV at the bottom left containing compound solid interfaces
$S_{\mathrm{cv,top}}$ |  circumferential surface of CV at the top containing interface type $12$
$S_{\mathrm{cv,fluid}}$ |  surface of CV within fluid phases
$t$ |  time
$\widehat{t}_{k}$ |  start of dissipation period $k$
$\widetilde{t}_{k}$ |  start of equilibrium period $k$
${\widehat{\Delta t}}_{k}$ |  duration of dissipation period $k$
$\widetilde{\Delta t}_{k}$ |  duration of equilibrium period $k$
$\overline{T}_{i}$ |  energy term $i$ corresponding to total analysis period
$\widehat{T}_{i,k}$ |  energy term $i$ corresponding to dissipation period $k$
$\widetilde{T}_{i,k}$ |  energy term $i$ corresponding to equilibrium period $k$
$V$ |  volume
$V_{\mathrm{cv}}$ |  volume of CV
$V_{\mathrm{cv,fluid}}$ |  volume of CV containing fluid phases
$\overleftarrow{V}$ |  volume swept by $S_{\mathrm{cv,bl}}$ over advance duration
$\overrightarrow{V}$ |  volume swept by $S_{\mathrm{cv,br}}$ over advance duration
$\boldsymbol{v}_{\mathrm{}}$ |  velocity
$v_{\mathrm{cap}}$ |  speed of fluid movement during dissipation events
$v_{\mathrm{cv}}$ |  speed of control volume advance
$\boldsymbol{v}_{\mathrm{cv}}$ |  velocity of control volume advance
$\widetilde{\boldsymbol{v}_{\mathrm{}}}$ |  velocity of fluids during equilibrium periods
$\widehat{\boldsymbol{v}_{\mathrm{}}}$ |  velocity of fluids during dissipation periods
$\widetilde{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}$ |  velocity gradient within fluids during equilibrium periods
$\widehat{\boldsymbol{\nabla}_{\mathrm{}}\boldsymbol{v}_{\mathrm{}}}$ |  velocity gradient within fluids during dissipation periods
$X_{\mathrm{cv}}$ |  distance travelled by CV during total analysis period
$x_{\mathrm{cv}}$ |  distance from start of CV travel
$\boldsymbol{x}_{\mathrm{cap}}$ |  location of dissipation event