# Contact Angle Hysteresis on Rough Surfaces Part II: Energy Dissipation via
Microscale Interface Dynamics

Pawan Kumar1  Dalton J. E. Harvie1  1Multiphysics Fluid Dynamics Group,
Department of Chemical Engineering, University of Melbourne, Parkville, VIC,
3010, Australia

###### Abstract

The wetting behaviour of surfaces is important for various applications like
super-hydrophobic surfaces, enhanced oil recovery, mining of metal ores and
anti-icing surfaces etc. For rough surfaces, which are the rule rather than
the exception, designing textured surfaces that have wetting properties
tailored to suit these applications generally involves either bio-mimicry or
trial and error. Existing wetting theories such as the well-known Wenzel
(1936) and Cassie & Baxter (1944) models do not predict wetting regimes and
importantly, only give a single equilibrium angle rather than a range of
stable contact angles (hysteresis) as observed in reality. In this paper, we
use a roughness-scale mechanical energy balance (derived in part I) combined
with simulations of micro-scale interface dynamics based on open-source
software (Surface Evolver (Brakke, 1992)) to calculate the energy dissipation
during the motion of an interface over a chemically homogeneous rough surface.
This dissipation is used to predict contact angle hysteresis (CAH) from
knowledge of just the surface roughness topography and equilibrium contact
angle. We simulate interface dynamics over a surface decorated with a periodic
array of round-edge square pillars and show that the energy dissipated varies
approximately as $\phi\ln\phi$ with the area fraction ($\phi$), and becomes
zero as $\phi\rightarrow 0$. The CAH predicted by our method is in good
agreement with the experimental results of Forsberg et al. (2010), especially
at low area fractions. We also compute CAH for an interface moving at 45° to
the surface periodicity direction, showing that the experimental measurements
are bracketed by the 0° and 45° advance direction results.

##  1 Introduction

The presence of roughness can significantly affect the wettability of a solid
(Butt et al., 2022). The water repelling and self-cleaning property of a lotus
leaf (Barthlott & Neinhuis, 1997), the ability of rice leaves to direct water
droplets in a specific direction (Zhu et al., 2010), water capture from the
air by a desert beetle (Parker & Lawrence, 2001) and the ability of a water
strider to walk on water (Feng & Jiang, 2006) are some of the examples where
this relationship has been used in the nature. In many engineering and
industrial applications this relationship plays a very important role as well.
The efficiency of industrial floatation cells (Chau et al., 2009) and
electronic printed devices (Meyers et al., 2008), passive oil-water separation
(Li et al., 2013), self-cleaning surfaces (Xu et al., 2016), stay dry fabrics
(Bae et al., 2009) and novel lab-on-a-chip devices (Jia et al., 2019) are some
of the examples where wetting is industrially important. The relationship
between roughness and wettability of a surface is, however, not fully
understood and we still generally rely on either trial and error or bio-
mimicking for designing functional wetting surfaces (Yan et al., 2011; Sun et
al., 2005; Zhai et al., 2006; Feng et al., 2007).

The spreading of a liquid on a surface in the presence of another immiscible
fluid is characterised by the angle that the interface between the two fluids
makes with the solid (measured from inside the fluid for which the contact
angle is being measured). This angle is called the contact angle and for ideal
surfaces (perfectly smooth, chemically homogeneous, rigid and inert to both
the fluids), is given by Young’s equation (Young, 1805)

$\cos\theta_{\rm{e}}=\left(\frac{\sigma_{\rm{2S}}-\sigma_{\rm{1S}}}{\sigma_{12}}\right).$
(1)

Here, $\theta_{\rm{e}}$ is the Young’s angle, and $\sigma_{\rm{1S}}$,
$\sigma_{\rm{2S}}$ and $\sigma_{12}$ are, respectively, the fluid-1/solid,
fluid-2/solid and fluid-1/fluid-2 interfacial tensions. The contact angle can
be measured experimentally and is usually measured at a small distance away
from the solid surface. We shall refer to this as the macroscopic contact
angle ($\theta_{\rm{m}}$). For ideal surfaces, the macroscopic contact angle
is the same as Young’s angle, however, this is not true for real surfaces. On
a real surface, $\theta_{\rm{m}}$ is measured between the fluid/fluid
interface and an average flat surface, which represents an apparent contact
angle (Marmur, 2006). The presence of chemical and/or physical heterogeneities
on real surfaces can cause $\theta_{\rm{m}}$ to be significantly different
from $\theta_{\rm{e}}$.

On a rough surface, $\theta_{\rm{m}}$ is observed to depend upon the wetting
history of the surface, and has different values depending upon whether a
particular state of the system (for example a particular volume of a droplet)
is obtained by advancing or receding the three-phase contact line (TPCL, where
the fluids 1, 2 and the solid surface meet). When the interface is made to
advance slowly111In the case of a droplet, either by increasing the volume or
by tilting the surface on which the drop is deposited. Or by lowering the
surface in the pool of liquid in the case of a Wilhelmy plate (Law & Zhao,
2016) type setting., the apparent contact angle increases until it reaches the
maximum value, at which the TPCL starts advancing. This angle is called the
advancing contact angle ($\theta_{\rm{a}}$). Similarly, when the interface is
made to recede, the TPCL starts receding only when the apparent contact angle
attains a minimum value, known as receding contact angle ($\theta_{\rm{r}}$).
The difference between the advancing and receding contact angles is known as
contact angle hysteresis
($\Delta\theta_{\rm{cah}}=\theta_{\rm{a}}-\theta_{\rm{r}}$). Contact angle
hysteresis (CAH) is an important parameter for characterising wetting on rough
surfaces, however, it is still not fully understood.

The earliest relationship between surface roughness and the macroscopic
contact angle ($\theta_{\rm{m}}$) was given by Wenzel (1936), as

$\cos\theta_{\rm{m}}=r\cos\theta_{\rm{e}}.$ (2)

Here, $r$ is the roughness, defined as the ratio of the actual to the apparent
surface area of the solid. The Wenzel equation is phenomenological in nature
based on the observation that the total surface area of a solid increases due
to roughness, and is valid only for the homogeneous (or Wenzel) wetting state
(i.e. the droplet is in complete contact with the solid surface and all the
surface crevices are filled without any trapping of the surrounding immiscible
fluid). There is another possible wetting state, known as the composite (or
Cassie or Fakir) wetting state, in which the surrounding fluid is trapped
inside the tops of the solid crevices. Here, $\theta_{\rm{m}}$ is governed by
Cassie-Baxter equation (Cassie & Baxter, 1944). Both Cassie-Baxter and
Wenzel’s equation predict a single equilibrium contact angle value, however,
because of $\rm{CAH}$, a rough surface exhibits a number of equilibrium
contact angles (Huh & Mason, 1977; Oliver et al., 1977; Ramiasa et al., 2014).
As such, these equations cannot predict the hysteresis on a rough surface
(Good, 1952; Kusumaatmaja & Yeomans, 2007; Pomeau & Vannimenus, 1985; Gao &
McCarthy, 2007; Patankar, 2010), and so are not very useful for designing
surfaces that require specific wetting behaviour.

Even though we lack a fully predictive wetting model for hysteresis, there are
a few fundamental works which can help in explaining the existence of
hysteresis on a rough surface. Shuttleworth & Bailey (1948) showed that a
droplet may exist in a number of metastable states on a rough surface, which
they suggested is responsible for the multiple equilibrium contact angles that
characterise CAH. Johnson and Dettre conducted numerical (Johnson & Dettre,
1964) and experimental (Dettre & Johnson, 1964) studies based on the idea of
the existence of multiple metastable states (Shuttleworth & Bailey, 1948).
They hypothesised that the origin of CAH lies in an energy barrier that must
be overcome to move between different metastable states and available droplet
vibrational energy is responsible for the contact angle hysteresis, however,
they did not present any method to quantitatively estimate the droplet
vibrational energy or relevant energy barriers. Cox (1983) developed a model
for hysteresis by treating the interface near the TPCL as a minimal surface
(zero mean curvature). The model was developed for sinusoidally undulating
surfaces that have a very small slope, which limits its application to real
surfaces which usually have more step-like topologies (Brandon et al., 2003;
Dettre & Johnson, 1964). Joanny & de Gennes (1984) presented a model for
contact angle hysteresis on a surface containing a dilute number of ‘strong
defects’ in the range of a few microns. They suggested that the TPCL pins on a
defect when the pinning force due to the defect and the elastic force acting
on the distorted TPCL balances out. There are multiple solutions to this force
balance and the TPCL can jump from one such position to the other, dissipating
energy in the process. This dissipation is related to the contact angle
hysteresis as

$\cos\theta_{\rm{r}}-\cos\theta_{\rm{a}}=nW_{\rm{d}},$ (3)

where $n$ is the density of the surface defects and ${W_{\rm{d}}}$ is the
energy dissipation during one complete hysteresis cycle around a single
defect. Many experimental studies have been done to quantitatively estimate
the energy dissipation during hysteresis (Priest et al., 2009; Patankar, 2010;
Ramos et al., 2003; Ramos & Tanguy, 2006; Reyssat & Quéré, 2009; Delmas et
al., 2011) however, to our knowledge there is no analytical or numerical model
available for calculating energy dissipation on generally structured surfaces.

Surface Evolver (SE) (Brakke, 1992) is open-source software that has been
extensively used for studying contact hysteresis on heterogeneous surfaces
(Dorrer & Rühe, 2008, 2007b; Forsberg et al., 2010; Dorrer & Rühe, 2007a;
Brandon et al., 1997; Chen et al., 2005). Forsberg et al. (2010) used SE to
calculate advancing contact angles during wetting on a surface with micron-
sized pillars of square cross-section. They simulated the liquid-air interface
between pillars and a flat surface parallel to the bottom of the pillars. For
obtaining the advancing contact angle, the apparent contact angle of the
interface was increased until it became unstable, either by depinning from the
pillars or by touching the neighbouring pillars. Semprebon et al. (2012) used
SE to simulate interface profiles pinned on cylindrical pillars and calculated
the advancing contact angle of the interface as the maximum possible apparent
angle for which an equilibrium interface morphology was possible, given the
constraints of equilibrium angle ($\theta_{\rm{e}}$), pillar geometry and
pillar separation distance. In all of these works (Dorrer & Rühe, 2008;
Forsberg et al., 2010; Semprebon et al., 2012) the roughness was attributed to
pillar structures (square or circular cross-section) with sharp edges, whereas
real surfaces comprise of defects with rounded edges. Promraksa & Chen (2012)
used SE to simulate the static equilibrium of a full droplet on a double
cosine wave-type surface. The advancing and receding contact angles here were
calculated as the maximum and minimum contact angles corresponding to certain
metastable states. Other than SE, Mognetti & Yeomans (2010) used Lattice
Boltzmann simulations to obtain the depinning angles for receding contact
lines on a superhydrophobic surface. In all these works, a static interface
shape or a static equilibrium shape of a full droplet was calculated rather
than the dynamics of the interface over a surface with ‘strong-defects’
(Joanny & de Gennes, 1984).

In the accompanying paper (part I) we developed a framework for predicting CAH
based upon a mechanical energy balance applied to a control volume moving with
the advancing TPCL. Key to this theory is being able to predict interface
dynamics around the TPCL from knowledge of the solid topology. In this work,
we use a novel numerical method to simulate the microscale dynamics of an
interface between two immiscible fluids moving at a very small velocity over a
rough solid. The numerical method is integrated into the mechanical energy
balance framework presented in part I to develop a fully predictive equation
capable of predicting contact angle hysteresis from just the knowledge of
surface topology and the equilibrium contact angle. The surface topology
considered in this work consists of periodically structured pillars of square
cross-section.

##  2 Theory

###  2.1 Dynamics of an interface advancing over a rough solid

In mechanical equilibrium, the Young-Laplace equation gives an expression for
the pressure difference ($\Delta p$) across an interface between two
immiscible fluids as

$\Delta p=2\kappa\sigma_{12},$ (4)

where $\kappa$ is the mean curvature of the interface, $\Delta\rho$ is the
density difference between the two fluids and $\sigma_{12}$ is the surface
tension of the fluid-fluid interface. For an interface between two immiscible
fluids moving over a rough surface ($h_{\rm{rough}}$) such that the length
scale of the surrounding fluid flow ($h_{\rm{surround}}$) is much greater than
the surface roughness ($h_{\rm{surround}}\gg h_{\rm{rough}}$), then in a
region close to the TPCL, $\Delta
p=O(\frac{\sigma_{12}}{h_{\rm{surround}}})\approx 0$. Putting $\Delta p=0$ in
equation (4) yields $\kappa=0$ which defines a minimal surface. Therefore, the
shape of the interface in equilibrium, near the TPCL can be represented by a
minimal surface (de Gennes et al., 2004, p. 13), provided that the length
scale of the surface roughness is very small as compared to that of the
surrounding fluid flow.

Turning to the dynamics of the interface, the augmented Navier-Stokes equation
for the flow of incompressible and immiscible fluids can be written as
(Popinet, 2009)

$\frac{D}{Dt}\rho\boldsymbol{v}=-\boldsymbol{\nabla}p+\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}+\sigma_{12}\kappa\delta_{\rm{S}}\boldsymbol{n}_{\rm{S}}+\rho\boldsymbol{g}.$
(5)

Here, $\boldsymbol{v}$ is the velocity, $\rho$ is the density, $p$ is the
pressure and $\mu$ is the dynamic viscosity of the fluid, $\sigma_{12}$ is the
surface tension of the interface, $\delta_{\rm{S}}$ is the surface delta
function, $\kappa=-\boldsymbol{\nabla}\cdot\boldsymbol{n}_{\rm{S}}$ is the
mean curvature of the interface, $\boldsymbol{n}_{\rm{S}}$ is the unit normal
to the interface directed into fluid-1 from fluid-2 and $\boldsymbol{g}$ is
the acceleration due to gravity. For the interface moving at very small
velocities ($\boldsymbol{v}\approx\boldsymbol{0}$) and in the absence of
pressure gradients ($\boldsymbol{\nabla}p=\boldsymbol{0}$), the Navier-Stokes
equation simplifies to $\kappa=0$ (upon neglecting the effect of gravity),
which is the same solution as given by the Young-Laplace equation. Indeed, we
say an interface is in equilibrium when two conditions are met: Namely that 1)
$\kappa=0$ everywhere on the interface, and 2) that the interface is locally
intersecting with the rough solid surface at Young’s angle. If either of the
conditions is not met then the interface will move, as we now explain.

If the interface is not in equilibrium due to the first condition not being
met, then considering the kinematics of the interface near the TPCL we will
find that $\kappa<0$ there. This is the first scenario shown in figure 2,
which illustrates how interface jumps occur. Violation of the first
equilibrium condition occurs when more liquid is added to a drop in
equilibrium or when a Wilhelmy plate is pushed into a pool of liquid, for
example. To understand the interface behaviour under these circumstances, we
refer back to the Navier-Stokes equation (5) and take the dot product of this
equation with the unit normal to the interface ($\boldsymbol{n}_{\rm{S}}$),
i.e.

$\left(\frac{D}{Dt}\rho\boldsymbol{v}\right)\cdot\boldsymbol{n}_{\rm{S}}=-\boldsymbol{\nabla}p\cdot\boldsymbol{n}_{\rm{S}}+\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}+\sigma_{12}\kappa\delta_{\rm{S}}\boldsymbol{n}_{\rm{S}}\cdot\boldsymbol{n}_{\rm{S}}+\rho\boldsymbol{g}\cdot\boldsymbol{n}_{\rm{S}}.$
(6)

Integrating equation (6) over a small control volume ($V_{\rm{CV}}$) around
the interface as shown in figure 1, and assuming for simplicity that
$\boldsymbol{n}_{\rm{S}}$ does not change significantly with time, gives

$\begin{split}\int_{{V_{\rm{CV}}}}\left(\frac{D}{Dt}\rho\boldsymbol{v}\right)\cdot\boldsymbol{n}_{\rm{S}}dV&=-\int_{{V_{\rm{CV}}}}\boldsymbol{\nabla}p\cdot\boldsymbol{n}_{\rm{S}}dV+\int_{{V_{\rm{CV}}}}\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}dV\\\
&+\int_{{V_{\rm{CV}}}}\sigma_{12}\kappa\delta_{\rm{S}}dV+\int_{{V_{\rm{CV}}}}\rho\boldsymbol{g}\cdot\boldsymbol{n}_{\rm{S}}dV.\end{split}$
(7)

The left-hand side of the equation (7) can be simplified using the Reynolds
transport theorem, noting that the control volume moves at the local velocity
($\boldsymbol{v}$), and assuming the fluid is incompressible,

$\int_{{V_{\rm{CV}}}}\frac{D}{Dt}\rho\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}dV=\frac{d}{dt}\int_{{V_{\rm{CV}}}}\rho(\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}})dV.$
(8)

Simplifying the pressure gradient term in equation (7) using the Gauss-
Ostrogradskii Divergence theorem to convert the volume integral to a surface
integral over the control volume surface ($S_{\rm{CV}}$), gives

$\int_{{V_{\rm{CV}}}}\boldsymbol{\nabla}p\cdot\boldsymbol{n}_{\rm{S}}dV=\int_{{S_{\rm{CV}}}}p(\boldsymbol{n}_{\rm{CV}}\cdot\boldsymbol{n}_{\rm{S}})dS=\Delta
pA,$ (9)

where $\Delta p$ is the pressure difference across the interface and $A$ is
the area of the interface. In equation (9) we have neglected pressure
contributions originating from the small (thickness $\epsilon$) sides of the
control volume. Using equations (8), (9) and (6) hence gives,

$\frac{d}{dt}\int_{{V_{\rm{CV}}}}\rho(\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}})dV=-\Delta
pA+\int_{{V_{\rm{CV}}}}\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}dV+\int_{{V_{\rm{CV}}}}\sigma_{12}\kappa\delta_{\rm{S}}dV+\int_{{V_{\rm{CV}}}}\rho\boldsymbol{g}\cdot\boldsymbol{n}_{\rm{S}}dV.$
(10)

We next assume for simplicity that the curvature of the interface in the small
$V_{\rm{CV}}$ is uniform, and denote the interfacial velocity in the direction
of $\boldsymbol{n}_{\rm{S}}$ as
$v_{\rm{S}}=\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}$ and the volume of the
control volume as depicted in figure 1 as $\delta V$. Also, since we are
interested in the interface dynamics near the TPCL, the system under
consideration is small (less than the capillary length $r_{\rm{cv,grav}}$ from
part I) and the effect of gravity ($\boldsymbol{g}$) can be neglected. This
gives

$\frac{d}{dt}(\rho v_{\rm{S}})\delta V=-\Delta
pA+\int_{{V_{\rm{CV}}}}\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}\cdot\boldsymbol{n}_{\rm{S}}dV+\sigma_{12}\kappa
A.$ (11)

This equation shows that the local acceleration of the interface is governed
by the pressure difference across the interface, the curvature of the
interface, and viscous dissipation occurring local to the interface. Noting
further based on the system geometry that the Laplace pressure difference
across the interface is given by $\Delta
p=\mathcal{O}(\sigma_{12}/h_{\rm{surround}})\approx 0$, equation (11) can be
further simplified to

$\frac{d}{dt}(\rho v_{\rm{S}})\delta
V=\int_{V_{\rm{CV}}}\left[\mu\boldsymbol{\nabla}^{2}\boldsymbol{v}\right]\cdot\boldsymbol{n}_{\rm{S}}dV+\sigma_{12}\kappa
A.$ (12)

This shows that for a negative $\kappa$, the acceleration of the interface in
the direction of $\boldsymbol{n}_{\rm{S}}$ is negative, i.e. the interface
accelerates in a direction opposite to the surface normal pointing from
fluid-2 to fluid-1. Likewise, for a positive $\kappa$, the interface
accelerates in the direction of $\boldsymbol{n}_{\rm{S}}$. The former scenario
represents an advancing interface while the latter refers to the receding
motion of the interface. Hence, equation (12) shows that if an interface is in
equilibrium ($\boldsymbol{v}=\boldsymbol{0}$) but is then perturbed to be just
out of equilibrium (so that $\kappa\neq 0$), the interface will accelerate in
the direction of $\kappa\boldsymbol{n}_{\rm{S}}$, and further the interface
will keep moving until a new equilibrium topology is found (at which
$\kappa=0$ again). More specifically, if the TPCL is pinned while the
macroscopic interface continues to advance, eventually $\kappa<0$ near the
TPCL and locally the interface will accelerate towards fluid 2 again to
restore equilibrium conditions.

Figure 1: Control volume of small thickness ($\epsilon\gg 1$) chosen around
the interface for integrating the Navier-Stokes equation. Figure 2:
Representation of the phenomena of TPCL jump when the interface is pushed
beyond a critical state. Scenario 1 shows the non-zero mean curvature of the
interface and scenario 2 represents the local contact angle
($\theta_{\rm{local}}$) being greater than Young’s angle ($\theta_{\rm{e}}$).
The outcome of either of the two scenarios is the jumping of the TPCL towards
the bulk of fluid 2, followed by subsequent pinning at the next defect in the
direction of the interface motion.

Under the second scenario the interface is not in equilibrium because the
local contact angle is not everywhere equal to the equilibrium angle, as
depicted by the second case in figure 2. Specifically, we assume that after
advancing the interface we have $\theta_{\rm{local}}>\theta_{\rm{e}}$, where
$\theta_{\rm{local}}$ is the local contact angle. This is another possibility
when more liquid is added to a drop which is in equilibrium or when a Wilhelmy
plate is pushed into a pool of liquid. Since $\theta_{\rm{local}}$ is not the
same as the equilibrium angle ($\theta_{\rm{e}}$), a driving capillary force
equalling $\sigma_{12}(\cos\theta_{\rm{e}}-\cos\theta_{\rm{local}})$ (Ramiasa
et al., 2014) accelerates the TPCL in the direction of interface advancement.

In summary, two conditions must be satisfied for an interface to be in
equilibrium (Bartell & Shepard, 1953), being 1) zero mean curvature
($\kappa=0$) everywhere on the interface and 2) the interface intersecting the
solid surface locally at Young’s angle
($\theta_{\rm{local}}=\theta_{\rm{e}}$). The above analysis shows however that
if either condition is not met then the interface will move at capillary
driven speeds in the direction of the macroscopic interface motion until a new
equilibrium position is found. We refer to this rapid motion of the interface
as a jump, and the overall behaviour of the TPCL as stick-slip. We denote the
state of an interface in equilibrium but where an infinitesimal movement of
the macroscopic interface will result in an acceleration of the interface as
the first critical state (that is, before the jump), and the first equilibrium
interface position found after this while the macroscopic interface continues
to move as the second critical state (that is, after the jump). Therefore, an
interface moving at infinitesimal small velocity under the action of an
infinitesimal pressure gradient can be modeled as a series of minimal surfaces
with occasional contact line jumps occurring between our defined first and
second critical states.

###  2.2 Simulating interfacial dynamics

With a physical model for interface advance established, we now consider how
these dynamics can be simulated.

####  2.2.1 Simulation domain

To simulate the interfacial dynamics we choose a simulation domain in the form
of a rectangular channel, as shown in figure 3(a). The bottom surface of the
simulation domain is decorated with pillars of square cross-section (side $a$
and height $h$) arranged in a square array with a spacing $d$ (see figure
3(b)). The height, width and length of the simulation domain are $H$, $W$ and
$L$ respectively. All lengths are non-dimensionalized with respect to the
pillar side, i.e.

$\begin{split}&W=\frac{W^{{}^{\prime}}}{a^{{}^{\prime}}},\quad
H=\frac{H^{{}^{\prime}}}{a^{{}^{\prime}}},\quad
L=\frac{L^{{}^{\prime}}}{a^{{}^{\prime}}},\\\
&d=\frac{d^{{}^{\prime}}}{a^{{}^{\prime}}},\quad
h=\frac{h^{{}^{\prime}}}{a^{{}^{\prime}}},\quad\text{and}\quad
a=\frac{a^{{}^{\prime}}}{a^{{}^{\prime}}}=1.\\\ \end{split}$ (13)

Here, $W^{{}^{\prime}}$, $H^{{}^{\prime}}$ and $L^{{}^{\prime}}$ are the
dimensional channel width, height and length respectively and
$a^{{}^{\prime}}$, $h^{{}^{\prime}}$, $d^{{}^{\prime}}$ are the dimensional
pillar side, height and inter-pillar distance respectively.

Figure 3: (a) Three-dimensional view of the simulation domain. Fluid-1 is
advancing towards the right while fluid-2 is receding. The macroscopic contact
angle ($\theta_{\rm{m}}$) is measured at the top of the domain. The symmetry
plane passing through the pillar centers and parallel to the $x-z$ plane is
shown in grey. (b) Advancing TPCL on a surface with roughness in the form of
micro-scale pillars (side $a$ and height $h$) arranged in a square array with
inter-pillar spacing $d$. The blue line shows the TPCL pinned at the pillars
at an instant of time ($t=\tau_{1}$) and the red line shows the TPCL pinned on
the next row of pillars at an instant of time ($t=\tau_{2}$) depicting the
TPCL jump. (c) Zoomed view of the interface near the TPCL capturing the
principal curvatures $\kappa_{1}$ and $\kappa_{2}$ in the $x-z$ and $x-y$
planes respectively.

Since the surface is periodic and the interface is moving in the direction of
surface periodicity, the equilibrium interface morphology will be periodic in
nature with the period the same as the inter-pillar distance ($d$). Why these
interface morphologies are symmetric is discussed later. For now, because of
this periodic nature, the width of the simulation domain ($W$) is set equal to
the distance between pillar centers ($d$) and the length of the simulation
domain is chosen to be long enough so as to capture at least one TPCL jumping
event after the simulation start. The height of the simulation domain is
chosen such that all the deformations in the interface are contained within it
(see sensitivity analysis in Appendix F). A typical simulation domain is
chosen such that the pillar is located at the center of the domain, i.e. at
$y=0$ and the domain walls are at $y=-W/2$ and $y=W/2$ respectively.

####  2.2.2 Energy minimization

As discussed in §2.1, physically the interface moves in one direction,
progressing through a series of equilibrium stages. These equilibrium
interface morphologies are obtained by minimising the total interfacial energy
within the simulation domain which comprises the interfacial energies of the
fluid-1/fluid-2, fluid-1/solid and fluid-2/solid interfaces respectively. The
total interfacial energy ($E^{{}^{\prime}}$) within the simulation domain can
be written as

$E^{{}^{\prime}}=\sum_{i<j}\sigma_{ij}A_{ij}=\sigma_{12}A_{12}+\sigma_{\rm{1S}}A_{\rm{1S}}+\sigma_{\rm{2S}}A_{\rm{2S}}.$
(14)

Here, $\sigma_{ij}$ and $A_{ij}$ represent the interfacial tension and area of
the $ij{\rm{th}}$ interface i.e., $12$, $\rm{1S}$ and $\rm{2S}$ representing
fluid-1/fluid-2, fluid-1/solid and fluid-2/solid interfaces. If the total area
of the bottom surface is $A_{\rm{t}}$, then
$A_{\rm{1S}}+A_{\rm{2S}}=A_{\rm{t}}$, i.e.
$A_{\rm{2S}}=A_{\rm{t}}-A_{\rm{1S}}$, and therefore equation (14) can also be
written as

$E^{{}^{\prime}}=\sigma_{12}A_{12}+(\sigma_{\rm{1S}}-\sigma_{\rm{2S}})A_{\rm{1S}}+\sigma_{\rm{2S}}A_{\rm{t}}.$
(15)

Defining $E=E^{{}^{\prime}}-\sigma_{\rm{2S}}A_{\rm{t}}$, where
$\sigma_{\rm{2S}}A_{\rm{t}}$ is a constant for a given fluid-2/solid
combination and simulation domain geometry, the total relative energy ($E$)
within the domain can be written using equation (15) and Young’s equation (1)
as

$E=\sigma_{12}(A_{\rm{12}}-A_{\rm{1S}}\cos\theta_{\rm{e}}).$ (16)

For the purpose of brevity in the rest of the paper we refer to the total
relative energy (equation (16)) as total energy.

The minimal surface morphology of interface is obtained by minimising the
total energy (equation (16)) within the simulation domain (as shown in figure
3(a)). The energy minimization is carried out in Surface Evolver (SE) (Brakke,
1992), which is an open-source software often used for studying surface
wetting problems (Irannezhad et al., 2023; Pour & Thiessen, 2019; Semprebon et
al., 2012; Forsberg et al., 2010; Dorrer & Rühe, 2008, 2007b, 2007a; Chen et
al., 2005; Brandon et al., 1997). SE uses the method of gradient descent for
arriving at the state of minimum energy. A surface is represented by a set of
triangular facets and the minimum energy configuration is obtained by moving
vertices in the direction of a negative gradient of energy while maintaining
the applied constraints of constant volume, fixed pressure, fixed edges,
vertices or any other user-defined constraints. An initial geometry is loaded
into the software, which then evolves towards a minimum energy shape with
subsequent refinement and iterations.

Figure 4: Use of superquadrics for modeling pillars. The shape of the cross-
section is fixed by the coefficient $\epsilon_{2}$, which is chosen as 4 for
square pillars. The sharpness of the edges can be controlled by the
coefficients $\epsilon_{1}$ and $\epsilon_{3}$. The higher the value of
$\epsilon_{1}$, $\epsilon_{3}$, the sharper the pillar edge. In this work we
use $\epsilon_{1}$=$\epsilon_{3}$=40 and $\epsilon_{2}$=4.

The SE algorithm minimises the total energy of the simulation domain as given
in equation (16), by adding the energy of individual facets of the interface
mesh via

$E=\left(\sigma_{12}\sum_{j=1}^{j={N_{1}}}\frac{1}{2}|\boldsymbol{e_{0}}\times\boldsymbol{e_{1}}|_{j}\right)_{12}-\cos\theta_{\rm{e}}\left(\sigma_{12}\sum_{k=1}^{{k=N_{2}}}\frac{1}{2}|\boldsymbol{e_{0}}\times\boldsymbol{e_{1}}|_{k}\right)_{1\rm{S}}.$
(17)

Here, $\boldsymbol{e_{0}},\boldsymbol{e_{1}}$ represents the facet edge
vectors and ${N}_{1},{N}_{2}$ represents the total number of facets on the
fluid-1/fluid-2 and fluid-1/solid interface respectively. The Surface Evolver
documentation (Brakke, 1994) gives more details regarding the implementation
of equation (17).

####  2.2.3 Pillar geometry

For modeling the pillars we use three-dimensional shapes known as
superquadrics (Barr, 1981): In our application, these are rectangular prisms
with rounded edges. In general, superquadrics can be represented by the
following equation

$\left(\frac{x-x_{1}}{s_{1}}\right)^{\epsilon_{1}}+\left(\frac{y-y_{1}}{s_{2}}\right)^{\epsilon_{2}}+\left(\frac{z-z_{1}}{s_{3}}\right)^{\epsilon_{3}}=1.$
(18)

where $\epsilon_{1},\epsilon_{2}$ and $\epsilon_{3}$ control the shape of the
superquadric cross-section, and $s_{1},s_{2},s_{3}$ control the aspect ratio
of the pillars respectively. The values of $x_{1},y_{1},z_{1}$ place the
pillar inside the domain. The main advantage of using superquadrics for
modeling pillar structures is that they create smooth surfaces with controlled
sharpness of the edges (see figure 4). As well as being numerically more
convenient (avoiding infinite curvatures associated with sharp corners), a
pillar with slightly rounded edges better represents the actual profile of
fabricated surfaces. For the current simulations we use
$\epsilon_{1}=40,\epsilon_{2}=4,\epsilon_{3}=40$, as shown in figure 4.

####  2.2.4 Simulating interface advancement

As already discussed in §2.1, an advancing interface moving at an
infinitesimal velocity under an infinitesimal pressure gradient can be
represented by a series of minimal surfaces, which are obtained by minimising
the total interfacial energy within the simulation domain (equation (16)) at
each simulated time point. Instead of the full interface, we simulate a
portion of a larger interface which is moving continuously under the action of
macroscopic flow at large distances away from the TPCL. From the point of view
of our simulations, this implies that the interface is moving at a constant
rate at the top of the domain (assuming the domain is sufficiently high such
that the separation of length scales inherent in the physical system is
ensured) whereas the TPCL moves in its characteristic stick-slip motion.
Therefore, the advancing motion of the interface within the simulation domain
can be attained by gradually advancing the interface top. In terms of boundary
conditions, this means that the interface top is constrained to a specific
value ($x_{\rm{top}}$) at a certain time point, which is then gradually
incremented in steps. The other boundary conditions are the orthogonality
between the simulation domain walls and the interface, i.e.

$\boldsymbol{n}_{\rm{S}}\cdot\boldsymbol{j}=0,$ (19)

along the domain walls, representing periodicity in the direction
perpendicular to the interface advancement direction, where
$\boldsymbol{n}_{\rm{S}}$ is the unit normal vector to the fluid-fluid
interface (directed into fluid-1 from fluid-2) and $\boldsymbol{j}$ is the
unit vector in the direction of $y$ axis (see figure 3(a)). Along the base of
the computational domain, the boundary condition of Young’s angle is applied,
that is, equation (1)

$\boldsymbol{n}_{\rm{S}}\cdot\boldsymbol{n}_{\rm{solid}}=-\cos\theta_{\rm{e}},$
(20)

where $\boldsymbol{n}_{\rm{solid}}$ is the unit normal vector to the solid
surface at the TPCL (directed out of the surface).

Figure 5 schematically illustrates our interface advance algorithm. A typical
simulation starts with a flat interface constrained according to the above-
mentioned boundary conditions and positioned close to a pillar such that the
interface is just touching it. We use a starting interface position
($x_{\rm{top}}$) such that the macroscopic contact angle ($\theta_{\rm{m}}$)
is equal to the Young’s angle ($\theta_{\rm{e}}$), and therefore, we choose
$x_{\rm{top}}=x_{0}-H/\tan\theta_{\rm{e}}$, where $x_{0}$ represents the TPCL
position. Now, the intersection between the interface mesh and the pillar is
identified. The intersecting facets of the interface and the associated edges
and vertices are constrained to follow the profile of the pillar (equation
(18)) and the surface energy of the intersecting facets (of the interface) are
changed to $-\sigma_{12}\cos\theta_{\rm{e}}$ (from $\sigma_{12}$) representing
their transition from a fluid-1/fluid-2 to fluid-1/solid interface. Now the
energy minimization is carried out in a number of steps with gradual mesh
refinements and an equilibrium interface morphology is obtained. This is the
initial equilibrium shape with $x_{\rm{top}}$ chosen such that
$\theta_{\rm{m}}=\theta_{\rm{e}}$. In order to advance the interface,
$x_{\rm{top}}$ is increased by a small amount ($\Delta x$), which is the step
size of the interface advancement. With $x_{\rm{top}}$ changed to
$x_{\rm{top}}+\Delta x$ at the top and keeping the other boundary conditions
unchanged, energy minimization is carried out once again to obtain the minimal
surface morphology of the interface at this new position. This process is
repeated until the interface reaches the first critical state, where it is not
possible to advance the interface and find a minimal surface shape such that
the TPCL is still pinned on the first pillar. The position of the interface
top and the macroscopic contact angle when the interface is in this critical
state are represented as $x_{\rm{s}}$ and $\theta_{\rm{s}}$ respectively. The
accuracy with which $x_{\rm{s}}$ and $\theta_{\rm{s}}$ can be predicted
depends upon the step size ($\Delta x$), with better accuracy achieved using
smaller step sizes (see Appendix A).

Figure 5: Flowchart of the algorithm used for simulating the advancement of an
interface by the incremental advance method.

Up until this point, we have discussed the advancing motion of an interface
starting from an equilibrium position with $\theta_{\rm{m}}=\theta_{\rm{e}}$
to the first critical state. Now, we discuss how the interface advance is
simulated after reaching this state. Any further increment in the
$x_{\rm{top}}$ now results in the depinning of the TPCL from the pillar222In
this paper, we focus on the interface motion such that the TPCL depins from a
pillar before it gets pinned again at the next suitable location and the
critical interface morphologies are such that the interface is not touching
the next pillar in the direction of the motion. and following our physical
model of interface advance the TPCL now moves forward until becoming pinned
again at the next pillar in the direction of advancement. To find an
equilibrium interface morphology, we start with finding the intersection
between the interface and the next pillar in the direction of advancement
followed by constraining the intersecting facets (and edges and vertices of
the interface) to the shape of the new pillar. The surface energy of the
intersecting facets of the interface is once again changed to
$-\sigma_{12}\cos\theta_{\rm{e}}$. With the interface touching the new pillar
and the above-mentioned changes to the intersecting facets of the interface,
energy minimization is carried out to obtain the equilibrium interface
morphology. We call this equilibrium interface morphology after the TPCL has
just depinned from a pillar as the ‘second-critical state’ and refer to the
corresponding interface position and the macroscopic contact angle as
$x_{\rm{s}^{{}^{\prime}}}$ and $\theta_{\rm{s}^{{}^{\prime}}}$ respectively.
The morphological transition between the two critical states captures a
typical TPCL jumping event. To further advance the interface, subsequent
$x_{\rm{top}}$ increments and energy minimization are carried out to obtain
the next equilibrium, first and second critical interface morphologies and
this process is repeated until we reach the desired maximum displacement of
the interface ($x_{\rm{max}}$). We call this method of simulating interface
dynamics by gradual increments in $x_{\rm{top}}$ the ‘incremental advance
method’. In our simulations, we were able to find an equilibrium interface
morphology on a neighbouring pillar in the direction of interface motion every
time the TPCL executes a jumping event as long as the step size is small
enough (see Appendix A). This is expected since the surface is structured and
all the pillars have the same geometry (Huh & Mason, 1977). Equilibrium
morphologies depicting a typical interface advancement are shown in Appendix
B.

Now that we have described the method of simulating a typical interface
advancement, we discuss the role of the step size ($\Delta x$) on the accuracy
with which the advancement of a real interface can be captured by our method.
A typical interface advancement on a rough surface is characterised by the
stick-slip motion executed by the TPCL. This has been discussed in detail in
part I and is briefly discussed here. The TPCL during a jumping event moves
with the capillary velocity ($\boldsymbol{v}_{\rm{cap}}$) and the jumping
event is completed in a very short duration of time ($\tau_{\rm{cap}}$), such
that

$v_{\rm{cap}}=\text{min}\left(\frac{\sigma_{12}}{\mu},\sqrt{\frac{\sigma_{12}}{\rho
d}}\right),$ (21)

and

$\tau_{\rm{cap}}=\frac{d}{v_{\rm{cap}}},$ (22)

where $\rho$ and $\mu$ are chosen as the maximum of the fluid’s densities and
viscosities (see part I). In a typical scenario of a water droplet spreading
in air on a surface where the roughness is of the order of a few microns,
$v_{\rm{cap}}$ is approximately 3 m/s and $\tau_{\rm{cap}}\approx$ 3 $\mu$s.
Therefore, during a TPCL jumping event, the surrounding interface doesn’t
significantly move macroscopically while the TPCL moves from one set of
pillars to the next. In the context of our simulations this means that during
the TPCL jump $x_{\rm{top}}$ should not change, i.e. both of the critical
states should be calculated using the same $x_{\rm{top}}$ (i.e.
$x_{\rm{s}^{{}^{\prime}}}\to x_{\rm{s}}$). In order to ensure this, we ideally
require an infinitesimally small step size to be used in our simulations,
(i.e. $\Delta x\to 0$), however this would result in a correspondingly large
number of equilibrium interface morphologies between the initial and the first
critical state and hence an impractically large computational time. However,
as we discuss in §2.3, for calculating the dissipation in energy during the
interface advancement, only the two critical states of the interface during
the advancing motion need be computed. Therefore, we propose another method
which is more efficient in finding these two states compared to the method of
incremental advance as described above. We call this method the ‘critical-
state method’, and it is described in Appendix C. In this method, the
interface top is constrained by a certain macroscopic angle
($\theta_{\rm{m}}$) and is advanced by gradually incrementing
$\theta_{\rm{m}}$, i.e. $\theta_{\rm{m}}+\Delta\theta$, where $\Delta\theta$
is the step size. For the case of an interface advancing on a structured
surface in the direction of surface periodicity, the critical-state method is
consistent with the incremental advance method (see Appendix D). Most of the
results presented in this paper are simulated using this critical-state
method.

In figure 6 we show example results for the typical advancement of an
interface between two immiscible fluids (fluid-1 and fluid-2) over a rough
surface similar to the one as shown in figure 3(b). We simulate the motion of
the interface in the direction of surface periodicity. Suppose at a particular
instant of time ($t=\tau_{1}$), the TPCL is shown by the solid blue lines in
figure 3(b). The interface is then made to advance continuously, such that at
the time ($t=\tau_{2}$), the TPCL is represented by the solid red line in
figure 3(b), after having executed a jump. Figure 6(a) shows the variation in
macroscopic contact angle ($\theta_{\rm{m}}$) with the interface position
($x_{\rm{top}}$) as it advances within the simulation domain.

Figure 6: (a) Variation in macroscopic contact angle ($\theta_{\rm{m}}$)
measured at the simulation domain’s top with the position of the interface top
($x_{\rm{top}}$) during advancement of the interface (for
$\theta_{\rm{e}}=72$°, $\phi=0.08$, $h/a=0.35$). Points ‘A’ and ‘B’ represent
the first and second critical states respectively. (b) TPCL lying on the
domain base corresponding to the equilibrium states shown in (a). (c)
Projection on the symmetry plane of the equilibrium interface morphologies
corresponding to the states shown in (a). The TPCL and the interface
projection on the symmetry plane corresponding to the first and second
critical states are shown in red in (b) and (c) respectively.

The portion of the TPCL lying on the domain base and the projection of the
interface on the symmetry plane (see figure 3(a)) for the interface
morphologies during the advancement of the interface is shown in figures 6(b)
and (c) respectively. For every value of $x_{\rm{top}}$ we have a unique
contact line shape. This is because as we change $x_{\rm{top}}$, the curvature
in $x-z$ plane, i.e. $\kappa_{1}$ (figure 3(c)) changes by a small amount.
This change in $\kappa_{1}$ is balanced by the same change in $\kappa_{2}$
(curvature in $x-y$ plane, see figure 3(c)), but in the opposite direction so
that the mean curvature remains zero. Here, both the curvatures, i.e.,
$\kappa_{1}$ and $\kappa_{2}$ are measured in the same neighbourhood,
therefore, everywhere on the interface, $\kappa_{2}$ adjusts so as to balance
the changes in $\kappa_{1}$ when $x_{\rm{top}}$ changes. However, when the
interface reaches a critical state (for example, state ‘A’ in figure 6(a)), it
is not possible to find a minimal surface shape under the constraint of the
local angle being equal to Young’s angle everywhere on the TPCL, and the
interface accelerates in the direction of $x_{\rm{top}}$ increment. When the
interface is at first critical state, the macroscopic angle attains its
maximum value ($\theta_{\rm{s}}$). Upon a small advancement of the interface
from this first critical state, the interface starts to accelerate but pins
again at the next closest pillar in the direction of advancement.

Indeed, on a structured surface with pillars of the same geometry and chemical
properties, an interface moving in the direction of surface periodicity is
able to find an equilibrium morphology on the next closest pillar every time
the TPCL executes a jump, provided that the macroscopic flow velocities are
very small ($\boldsymbol{v}_{\rm{CV}}\approx\boldsymbol{0}$) and
$\Delta\theta$ is small. This is illustrated in figure 7, where we show the
equilibrium interface morphologies during the two critical states. The
interface motion during the TPCL jump behaves as if it is pivoted at the top
of the simulation domain and this motion results in a decrease in the
interface’s curvature in the $x-z$ plane. The interface remains pinned on this
new defect until the curvature in the $x-z$ plane is increased to the value at
which the interface arrives at the first critical state again and the TPCL
executes a jump again. This is repeated again and again as the interface
advances over the surface. The TPCL moves from one row of pillars to the next
row, pinning every time it encounters a fresh row. This motion captures in the
well-known stick-slip behaviour of the TPCL.

Figure 7: Equilibrium interface morphologies capturing a TPCL jump. The
macroscopic contact angle decreases from $\theta_{\rm{s}}$ to
$\theta_{\rm{s}^{{}^{\prime}}}$ during the morphological transition from the
critical to the second-critical state. The difference between the two angles
($\Delta\theta_{\rm{d}}=\theta_{\rm{s}}-\theta_{\rm{s}^{{}^{\prime}}}$)
decreases as the domain height is increased (see §3.3).

Another point regarding interface movement is that the morphology of an
equilibrium interface depends upon its history. Specifically, interface
behaviour during advancing and receding motion is different (hysteresis). Here
we show that the numerical model also predicts different interface
morphologies depending on its history. We start with the interface in
equilibrium and instead of advancing it, we recede the interface gradually (by
gradual decrements in $\theta_{\rm{m}}$). In figure 8(a) we plot the portion
of the TPCL which is on the domain base capturing the first and second
critical states during the advancing and receding motion of the interface. We
observe that these two critical states are completely different during the
advancing and receding motion. In figure 8(b) we show the projection of
interface morphologies on the symmetry plane capturing the first and second
critical states, which also shows the difference in the interface morphologies
during advancing and receding motion. Therefore, the interface morphologies
and the macroscopic contact angle in critical states not only depend upon the
surface roughness and Young’s angle but also upon the history of the
interface. This is consistent with our physical interpretation, discussed in
§2.1 regarding how an interface advances.

Figure 8: (a) TPCL on the domain base capturing first and second critical
state interface morphologies during advancing (red) and receding (blue)
motion. (b) Projection of the equilibrium interface morphologies on the
symmetry plane during advancing (red) and receding (blue) motion capturing
first and second critical states. The interface is traversing on a surface
with $\theta_{\rm{e}}=72$° and $h/a=0.35$.

###  2.3 Calculating energy dissipation from interface dynamics

Now that we have established an efficient numerical method for predicting
interfacial dynamics, we want to use it to find the dissipation in energy
($D$) as the interface advances within the simulation domain. Referring back
to part I, this dissipation in energy is directly related to CAH, and is given
by

$D=-\sum_{k=1}^{N}\sum_{i<j}\sigma_{ij}\frac{\widehat{\Delta
A_{ij}}_{k}}{A_{\rm{CV}}},$ (23)

where $\widehat{{\Delta A_{ij}}_{k}}$ is the change in the area of the $ij$
interface within the simulation domain when the TPCL is executing its $k$th
jumping event, $\sigma_{ij}$ is the interfacial tension of the $ij$ interface
and $A_{\rm{CV}}$ is the area traversed by the TPCL projected on the $x-y$
plane during which it executes $N$ number of jumping events.

Before going into the details of this energy dissipation calculation, we first
examine some aspects of the nature of total interfacial energy variation as an
interface advances within the simulation domain. Figure 9 shows the variation
in total non-dimensional energy ($\overline{E}$) within the simulation domain
as a function of $x_{\rm{top}}$. The total energy is non-dimensionalized by
the cross-sectional area of the pillars and the surface tension of the fluid-
fluid interface i.e.

$\overline{E}=\frac{E}{\sigma_{12}a^{2}}.$ (24)

As the interface advances within the domain, we observe a decreasing trend in
the total non-dimensional energy. This is because of the particular nature of
the surface under consideration that has $\theta_{\rm{e}}<90$°.

Figure 9: Variation of total non-dimensionalized energy within the simulation
domain ($\overline{E}$, shown by blue circles) with interface position
($x_{\rm{top}}$). The advancing motion of the interface has been captured on a
surface with $\theta_{\rm{e}}=72$°, $\phi=0.08$ and pillar aspect ratio of
0.35. As the interface advances, the TPCL moves in a stick-slip fashion
executing jumps and dissipating energy. Here we have shown three such
dissipation events. We observe that during the motion of an interface in the
surface periodicity direction of a structured surface, the magnitude of non-
dimensional energy dissipation ($\overline{D}_{1}$, per pillar) is the same
during each TPCL jumping event. Insets (A)-(C) show different equilibrium
morphologies of the interface during advancement. Inset (A) shows the
equilibrium interface morphology when $\theta_{\rm{m}}=\theta_{\rm{e}}$. This
represents a typical starting point of an interface advancement simulation.
Inset (B) shows the equilibrium interface morphology during a critical state
(the interface is advancing from the left to the right side). Inset (C) shows
the equilibrium interface morphology during a second-critical state. Insets
(B) and (C) represent the first TPCL jumping event. Variation in macroscopic
contact angle ($\theta_{\rm{m}}$, in degrees) corresponding to the equilibrium
states shown in the total non-dimensional energy ($\overline{E}$) plot is
shown by black squares. When the interface reaches a critical state (for
example state ‘B’), the macroscopic contact angle attains the maximum value
($\theta_{\rm{s}}$).

However, there is a sudden drop in the total non-dimensional energy at certain
points, such as ’B’ and ’C’. We also observe a drop in the macroscopic contact
angle at the same times. These points represent first and second critical
states respectively (see figure 6). Also in figure 9, we observe that the drop
in total non-dimensional energy (dissipation) is the same during each TPCL
jumping event. This is because the dynamics of each jump are the same as the
interface is moving over a structured surface and in the direction of the
surface periodicity.

Figure 10: (a) Three-dimensional view of the simulation domain with three
pillars in a row ($\phi=0.13,h/a=0.35,\theta_{\rm{e}}=72$°). The interface is
pinned on the first row and represents the stable morphology just before
depinning from the pillars (critical state). (b) Zoomed-up view of the pinned
interface. Equilibrium interface morphology is symmetric around each pillar.
The zoomed is flipped about the $z$ axis to show the fluid-1 tongue protruding
between the pillars.

We also examine the role of the simulation domain’s width in calculating
interface dynamics. Simulations can be performed in a computational domain of
any width which is an integer multiple of the inter-pillar distance. In figure
9 we have shown the interface dynamics calculated using a simulation domain
which has a width equal to the distance between neighbouring pillars ($W=d$),
while in figure 10 we show the equilibrium morphology of an interface pinned
on a row of three pillars. The shapes are the same. Further, figure 11 shows
the interface profile ($\phi=0.13$, $\theta_{\rm{e}}=72$° and $h/a=0.35$)
projected on a plane parallel to the $x-z$ plane and passing through $y=-W/2$,
$y=-W/4$ (i.e., passing through section BB in figure 10(b)), $y=W/4$ (i.e.,
passing through section AA in figure 10(b)) and $y=W/2$ from the three pillar
width simulation. All four profiles are identical, confirming the symmetry in
interface morphology around each pillar.

Figure 11: (a) Equilibrium interface morphologies projected on a plane
parallel to $x-z$ plane and passing through domain’s front and back walls and
through section ‘AA’ and ‘BB’ as shown in figure 10(b). Similar profiles are
obtained at the four locations. (b) Total non-dimensional energy dissipation
($\overline{D}A_{\rm{CV}}$) divided by the number of pillars along the row
($N_{y}$) ($h/a=0.35,\theta_{\rm{e}}=72$°). The energy dissipation per pillar
is independent of the number of pillars used in the simulation.

Another aspect of this periodicity in interface morphologies is that the
dissipation in energy from an interface pinned on a single pillar is the same
as the energy dissipation from an interface pinned on a row of $N_{y}$ pillars
divided by the number of pillars ($N_{y}$). In order to verify this, we
calculate the non-dimensional energy dissipation ($\overline{D}$, as the
difference in total non-dimensional interfacial energy within the domain
between the first and second critical states) when an interface pinned on a
row of three pillars along the simulation domain’s width executes a jump. The
dissipation is calculated for a single TPCL jumping event from one row to the
next (3 pillars in each row). We plot the total dissipation divided by the
number of pillars with the pillar area fraction ($\phi$) in figure 11(b) for
both the single pillar and three pillar width simualtions. We observe that for
a structured surface (and when the interface is moving in the direction of
surface periodicity) the dissipation per pillar is independent of the number
of pillars used in the simulation.

Hence, the energy dissipated when the TPCL pinned on a single pillar executes
a jump to the next pillar can be used to calculate the total dissipation in
energy during the travel of an interface that executes $N$ jumps over an array
of pillars that is $N_{y}$ pillars wide. Specially for this situation equation
(23) becomes

$D=-N_{y}N\sum_{i<j}\sigma_{ij}\frac{\widehat{\Delta
A_{ij}}_{1}}{A_{\rm{CV}}},$ (25)

where $\widehat{\Delta{A_{ij}}_{1}}$ is the change in interfacial areas during
a single TPCL jumping event from one pillar to the next. Equation (25) can,
therefore be written as

$\begin{split}D&=-N_{y}N\left(\frac{\sigma_{12}\Delta
A_{12}+\sigma_{\rm{1S}}\Delta A_{\rm{1S}}+\sigma_{\rm{2S}}\Delta
A_{\rm{2S}}}{A_{\rm{CV}}}\right),\\\ &=-N_{y}N\left(\frac{\sigma_{12}\Delta
A_{12}-(\sigma_{\rm{2S}}-\sigma_{\rm{1S}})\Delta
A_{\rm{1S}}}{A_{\rm{CV}}}\right).\end{split}$ (26)

In equation (26) we have used the relationship $\Delta A_{\rm{2S}}=-\Delta
A_{\rm{1S}}$. Using equations (16), (1) and (26), we can write

$\begin{split}D&=-N_{y}N\frac{\Delta E}{A_{\rm{CV}}}\\\
&=\sigma_{12}\phi\overline{D}_{1},\end{split}$ (27)

where $\phi=N_{y}NA_{\rm{p}}/A_{\rm{CV}}$ is the pillar area fraction,
$A_{\rm{p}}$ is the cross-sectional area of the pillar tops ($a^{2}$ in the
present case) and $\overline{D}_{1}$ is the non-dimensional dissipation per
pillar per dissipation event, defined by

$\overline{D}_{1}=\frac{-\Delta E}{\sigma_{12}A_{\rm{p}}}$ (28)

From equations (27) and (28) we can relate $D$, which is the total dissipation
in energy per unit area traversed by the TPCL ($A_{\rm{CV}}$) during interface
advancement over a structured surface in the direction of surface periodicity,
to the non-dimensional energy dissipation per pillar per dissipation event as
calculated in single pillar column simulations. The total dissipation
($\overline{D}$) can also be non-dimensionalized by the interfacial tension
$\sigma_{12}$, i.e.

$\overline{D}=\frac{D}{\sigma_{12}}=\phi\overline{D}_{1}.$ (29)

In the rest of the paper, we refer to the dissipation per pillar per
dissipation event as dissipation per pillar for brevity.

We finish our description of the numerical method with two final comments
regarding efficiency. As well as the periodicities discussed above, the fluid-
fluid interface also possesses a reflection symmetry with respect to a plane
passing through the centre of the pillars and parallel to the $x-z$ plane.
Mathematically this is expressed as
$\boldsymbol{n}_{\rm{S}}\cdot\boldsymbol{j}=0$ along this plane (figure 3(a),
symmetry plane is shown in grey colour). This symmetry allows us to reduce the
computational expense by simulating only half of the simulation domain width
(compared to a single pillar simulation) without affecting the outcome of the
simulations. Therefore, we use a domain which is only $d/2$ wide and $30d$
high for the majority of our simulations ((see Appendix F)). The dissipation
for the full domain (i.e., $H=30d,W=d$) is then obtained as twice the
dissipation calculated using this half domain (i.e., $H=30d,W=d/2$). Also, to
increase computational efficiency and accuracy the simulations are performed
using a non-uniform mesh which is dynamically adjusted based on distance to
the rear face of each pillar, according to equation (44). Appendix E gives
more details of this dynamic mesh generation technique.

##  3 Results and discussions

In this section, we present our results for the advancing and receding motion
of an interface on a surface with a structured array of square pillars. We
discuss the effect of pillar area fraction ($\phi$) and aspect ratio ($h/a$)
on the macroscopic contact angle and dissipation in energy. Based on the
numerical results, we develop a relationship between the pillar area fraction
and the energy dissipation for both the advancing and receding interface.
Finally, using the mechanical energy balance, we present a predictive equation
for the contact angle hysteresis on such surfaces.

###  3.1 Advancing interface

####  3.1.1 Interface behaviour

Figure 12(a) shows the equilibrium interface profile projected on the symmetry
plane (passing through pillar center and parallel to the $x-z$ plane) for a
pinned interface with $\theta_{\rm{e}}=72$°, $h/a=0.50$, $\phi=0.03,0.04,0.08$
and $0.13$. The solid lines represent the interface profile at the first
critical state and the dashed lines represent the equilibrium interface
profiles after the TPCL jump i.e. the second-critical state. We observe that
$\theta_{\rm{s}}$ increases with the increase in pillar area fraction and from
figure 12(b), where we have plotted the TPCL lying on the domain base, we can
see that the curvature in the $x-y$ plane adjacent to the pillar increases
with the increase in pillar area fraction ($\phi$).

Figure 12: (a) Equilibrium interface profiles projected on the symmetry plane
passing through the center of the pillar and parallel to $x-z$ plane. The
results are calculated for $\theta_{\rm{e}}=72$° and pillar aspect ratio
($h/a$) of 0.5. The critical contact angle ($\theta_{\rm{s}}$) exhibited by
the pinned interface varies with area fraction ($\phi$) ($110$°, $102$°,
$94$°, $89$° for pillar area fractions $\phi=0.13$ (red), $\phi=0.08$ (green),
$\phi=0.04$ (blue), $\phi=0.03$ (black) respectively). (b) Shows the TPCL on
the domain bottom corresponding to the equilibrium interface morphologies as
depicted in (a). The TPCL curvature increases as the pillar area fraction
increases.

Figure 13(a) shows the interface morphology during a first critical state
projected on the simulation domain wall ($y=W/2$) as well as the symmetry
plane for $\theta_{\rm{e}}=72$°, $h/a=0.35$ and for a large variety of area
fractions between $\phi=0.01$ to $\phi=0.70$. The maximum macroscopic contact
angle ($\theta_{\rm{s}}$) for different pillar area fractions ($\phi$) is also
shown. As previously noted, there is a direct relationship between
$\theta_{\rm{s}}$ and $\phi$, i.e. $\theta_{\rm{s}}$ increases with $\phi$.
Other information that can be inferred from figure 13(b), is the nature of the
projection on the domain wall. Since the interface in equilibrium represents a
minimal surface, at every point the mean curvature is zero. Therefore, a high
curvature of the interface in the $x-y$ plane is balanced by a similar
magnitude of curvature, but in the opposite direction, in the $x-z$ plane,
with both measured in the same neighbourhood. As noted earlier, the curvature
in the $x-y$ plane next to the pillar increases with $\phi$, therefore, the
curvature of the interface in the $x-z$ plane as well increases with $\phi$,
which is the reason we observe large $\theta_{\rm{s}}$ values at higher area
fractions. A high curvature in the $x-y$ plane also means that the TPCL
intersects the domain wall at a point closer to the pillar on which the
interface is pinned. Therefore, in figure 13(b) the interface projection on
the domain front wall for the pillar area fraction ($\phi=0.70$) intersects
the domain bottom very close to the pillar on which the interface is pinned
(this can also be observed in the TPCL projections for the first critical
interface morphologies shown as solid lines in figure 12(b)).

Figure 13: Equilibrium interface morphologies just before depinning (first
critical state) as projected on the symmetry plane (shown by dashed-dotted
lines) and the domain wall ($y=W/2$, shown by solid lines) for an advancing
interface on a surface with $\theta_{\rm{e}}=72$°, $h/a=0.35$ and area
fractions $\phi=0.01,0.03,0.04,0.08,0.13,0.20,0.30,0.40,0.50,0.70$. (b) Zoomed
view of interface morphologies near the TPCL.

Apart from the pillar area fraction, we observe that the pillar aspect ratio
also affects the interface curvature and the maximum macroscopic contact angle
($\theta_{\rm{s}}$). In figure 14 we show the projection of interface
morphologies in the first critical state for a pillar aspect ratio ($h/a$) of
1.5 ($\theta_{\rm{e}}=72$° and $\phi=0.01,0.03,0.04,0.08$ and $0.13$). We
observe that qualitatively, the interface in equilibrium behaves in a similar
manner on both the pillar aspect ratios ($h/a=0.35,1.5)$, however
quantitatively the principal curvatures in both the $x-z$ and $x-y$ planes are
comparatively higher for higher aspect ratio pillars.

Figure 14: Equilibrium interface morphologies just before depinning (first
critical state) as projected on the symmetry plane (shown by dashed-dotted
lines) and the domain wall ($y=W/2$, shown by solid lines) for an advancing
interface on a surface with $\theta_{\rm{e}}=72$°, $h/a=1.5$ and area
fractions $\phi=0.01,0.03,0.04,0.08,0.13$. (b) Zoomed view of interface
morphologies near the TPCL.

This is reflected in a higher magnitude of the maximum macroscopic contact
angle ($\theta_{\rm{s}}$) for higher aspect ratios which is shown in figure
15, where we plot the first critical interface morphologies projected on the
symmetry plane for pillar aspect ratios 0.35 and 1.5 and area fractions 0.01,
0.03, 0.04, 0.08 and 0.13 respectively.

Figure 15: First critical equilibrium interface morphologies as projected on
the symmetry plane for pillar aspect ratio ($h/a$) 0.35 (a) and 1.5 (b) for an
advancing interface on a surface with $\theta_{\rm{e}}=72$° and area fraction
$\phi=0.01,0.03,0.04,0.08,0.13$. The interface curvature as well as the
maximum macroscopic contact angle ($\theta_{\rm{s}}$) increases with the
pillar aspect ratio ($h/a$).

An interesting equilibrium interface morphology arises when the pillar area
fraction is increased beyond a certain critical point, at which we observe
that the interface makes contact with the next pillar before depinning from
the previous pillar (shown in figure 16(a)). The area fraction at which this
morphological transition occurs depends upon the pillar geometry and Young’s
angle. Even though such interface morphologies are not the focus of this
paper, we show one case as an example of the capability of the present
numerical method to deal with such morphologies. In figure 16 we show the
advancing motion of an interface over a surface with $\theta_{\rm{e}}=121$°
and $h/a=1.5$.

Figure 16: Equilibrium interface morphologies capturing the advancing motion
when the interface while pinned on a pillar touches the next pillar in the
direction of advancement. (a) The interface is touching the next pillar and
the fluid-1 spontaneously spreads on the top of the pillar as shown in (b).
However, the TPCL remains pinned on the first pillar while the fluid-1
gradually wets the pillar tops. During this process, a thin film of the
surrounding fluid (fluid-2) gets trapped in between the interface and the
solid surface as shown in (c). The interface morphologies shown here are
simulated with $\theta_{\rm{e}}=121$°, $\phi=0.13$ and $h/a=1.5$.

The interface, before depinning from the first pillar makes contact with the
next pillar in the direction of interface advancement and a portion of the
fluid-1/fluid-2 interface is converted to the fluid-1/solid interface. The
equilibrium interface morphology is shown in figure 16(b) where we observe
that fluid-1 only wets the pillar top. This process is repeated again,
resulting in an equilibrium interface morphology with fluid-1 wetting the
pillar tops while fluid-2 is trapped between the solid surface and the
interface (figure 16(c)). It is also possible to have other equilibrium
interface morphologies, such as the fluid-1 wetting the pillar tops as well as
some portion of the pillar sides for a different set of $\theta_{\rm{e}}$ and
$h/a$ values. Essentially this is modeling the transition from a Wenzel to
Cassie wetting state (for certain Young’s angles). The study of such cases is
left to future work.

####  3.1.2 Energy dissipation during advancing motion of the interface

In this section, we develop a relationship between the energy dissipation and
pillar area fraction for an interface advancing in the direction of surface
periodicity. Figure 17 shows the non-dimensional energy dissipation per pillar
($\overline{D}_{1}$) and total non-dimensional energy dissipation
($\overline{D}$) plotted against pillar area fraction ($\phi$) (pillar aspect
ratio 0.35 and 1.5, and $\theta_{\rm{e}}=72$°).

Figure 17: Variation in the total non-dimensional energy dissipation
($\overline{D}$, shown by triangles) and non-dimensional energy dissipation
per pillar ($\overline{D}_{1}$, shown by circles) with the pillar area
fraction ($\phi$) for an interface advancing in the surface periodicity
direction is shown. Results have been simulated for Young’s angle $72$° and
pillar aspect ratio of (a) 0.35 and (b) 1.5 respectively. Fitting equations of
dilute (equation (32)) and non-dilute (equation (34)) form for $\overline{D}$
are shown by dashed-dotted and solid black lines respectively. Fitting
equations for $\overline{D}_{1}$ are shown by dashed-dotted (dilute form,
equation (30)) and solid (non-dilute, equation (33)) blue lines respectively.

The non-dimensional dissipation per pillar is obtained from the interface
dynamics simulations for different area fractions and pillar aspect ratios,
while the total non-dimensional dissipation is obtained using,
$\overline{D}=\phi\overline{D}_{1}$ from equation (29).

For fitting $\overline{D}_{1}$ to an appropriate correlation, we start with an
equation form proposed by Joanny & de Gennes (1984) for the dissipation in
energy due to a single strong defect on a dilute surface, i.e.

$\overline{D}_{1}=A\ln\phi+C.$ (30)

Here $A$ and $C$ are constants that depend upon the pillar aspect ratio
($h/a$) and Young’s angle ($\theta_{\rm{e}}$) of the surface. The form of the
relationship between $\overline{D}_{1}$ and $\phi$ (equation (30)) arises from
the relationship between $\phi$ and $d$ (inter pillar distance) on a
structured surface (Joanny & de Gennes, 1984). When a liquid surface is
perturbed by the presence of a pillar, the perturbation is carried to
infinity, which can be captured by a function of logarithmic nature. Under
these conditions, the geometry of the pillar is of secondary importance to the
distance between pillars, provided that the pillar can actually pin the
interface and the distance between pillars is relatively large (i.e., a dilute
surface). The interface is perturbed when it touches a pillar, however, the
perturbation is arrested by the surrounding pillars which impose a cutoff
length on the deformation caused in the interface. Measured with respect to
the undistorted TPCL profile, if $x$ is the position of the distorted TPCL in
the direction of interface advance at a certain location ($y$) along the
domain width, then it is related to the inter-pillar distance ($d$) and the
location of the pillar center (i.e., $y_{1}$, see equation (18)) as (Nadkarni
& Garoff, 1992)

$x=\lambda\ln\left(\frac{d}{|y-y_{1}|}\right),$ (31)

where $\lambda$ is a constant depending upon the pillar geometry, equilibrium
angle ($\theta_{\rm{e}}$) and fluid/fluid interfacial tension ($\sigma_{12}$).
Also, since the pillar area fraction is a function of the inter-pillar
distance (that is, $\phi=1/d^{2}$, noting that $d$ is non-dimensionalised by
$a$), the TPCL profile of a pinned interface can be represented by the
logarithmic function of the pillar area fraction ($\phi$). This is shown in
figure 12(b), where we have plotted the TPCL for different area fractions. As
the energy dissipation depends upon the change in fluid-fluid and fluid-solid
interfacial areas, which varies logarithmically with $\phi$,
$\overline{D}_{1}$ displays a logarithmic dependence on $\phi$. Hence, under
the dilute surface assumption, the total non-dimensional dissipation
($\overline{D}$) can be obtained using equations (30) and (29) as

$\overline{D}=A\phi\ln\phi+C\phi$ (32)

At higher $\phi$, however, the relative distance between pillars is less and
the shape of the pillars becomes more important in determining the form of the
dissipation and we expect $\overline{D}$ to deviate from the dilute form.
Therefore we incorporate a term $B\phi$ into the expression, rationalised as
the first Taylor series correction to the dilute form as an alternative
correlation form for $\overline{D}_{1}$,

$\overline{D}_{1}=A\ln\phi+B\phi+C.$ (33)

Here, $B$ is a constant depending upon the pillar geometry and Young’s angle.
Based on equation (33), a ‘non-dilute’ form of the total non-dimensional
dissipation can similarly be written as

$\overline{D}=A\phi\ln\phi+B\phi^{2}+C\phi$ (34)

We calculate the coefficients $A$, $B$ and $C$ by fitting equations (32) and
(34) to the variation in total non-dimensional dissipation ($\overline{D}$)
with pillar area fractions generated from our numerical results. Here,
$\overline{D}$ is calculated from the interface dynamics simulations over a
single pillar as $\overline{D}=\phi\overline{D}_{1}$. The values of different
fitting coefficients and the $R^{2}$ values of the fit (for $\overline{D}$)
are given in table 1. Overall the dilute and non-dilute correlation forms fit
the results accurately, with the non-dilute form producing slightly higher
$R^{2}$ values.

Equation |  $h/a=0.35$ |  $h/a=1.5$
---|---|---
$\overline{D}=A\phi\ln\phi+C\phi$ (dilute surface) |  $A=-1.39$ $C=1.64$ $R^{2}=0.996$ |  $A=-7.44$ $C=-6.49$ $R^{2}=0.987$
$\overline{D}=A\phi\ln\phi+B\phi^{2}+C\phi$ (non-dilute surface) |  $A=-1.80$ $B=1.09$ $C=0.77$ $R^{2}=0.998$ |  $A=-11.26$ $B=53.16$ $C=-20.95$ $R^{2}=0.999$
Table 1: Fitting parameters for the total non-dimensional energy dissipation
($\overline{D}$) dependence on the pillar area fraction ($\phi$) for an
advancing interface ($\theta_{\rm{e}}=72$°), as suggested in equations (32)
and (34).

An observation from the peculiar nature of the $\overline{D}_{1}$ variation in
equations (30) and (33), is that despite the presence of $\ln\phi$, which
makes the dissipation per pillar non-vanishing even on very dilute surfaces
($\phi\to 0$), the behaviour of total dissipation ($\overline{D}$) is
different, becoming zero as the pillar area fraction approaches zero. In fact,
as shown by the simulation results $\overline{D}_{1}$ becomes very large as
the pillar area fraction is reduced, as can be seen in figure 17. The total
non-dimensional energy dissipation ($\overline{D}$) becomes zero on very
dilute surfaces due to the limit of $\phi\ln\phi$, i.e.

$\begin{split}\lim_{\phi\to 0}\overline{D}&=\lim_{\phi\to
0}(A\phi\ln\phi+B\phi^{2}+C\phi),\\\ &=\lim_{\phi\to 0}(A\phi\ln\phi),\\\
&=\lim_{\phi\to
0}\frac{A\ln\phi}{1/\phi}\stackrel{{\scriptstyle\text{H}}}{{=}}(A/\phi)(-\phi^{2}),\\\
&=0.\end{split}$ (35)

Based on the nature of $\overline{D}$ and $\overline{D}_{1}$ with respect to
$\phi$, it may appear that the behaviour of an interface on a surface with a
dilute number of strong defects is very different from the surface having a
single strong defect. According to equation (33), if an interface moving with
a small macroscopic flow velocity ($\boldsymbol{v}_{\rm{CV}}$) interacts with
a single strong defect then it should dissipate an infinite amount of energy
when released from the defect. However, in reality, physical surfaces have
finite sizes and (usually) more than a single defect, so have finite cutoff
lengths. For example for a small droplet, the cutoff length is the droplet
size. If a droplet of a few millimeters in diameter spreads on a surface with
a single strong defect (Joanny & de Gennes, 1984) of a few micrometers in
size, the area fraction would be of the order of 10-6. The non-dimensional
dissipation per pillar, according to the equation (33) and values as given in
table 1 is then of the order of 10-5, which is almost zero, as predicted by
equation (34). Also, from our simulations we observe that the maximum
macroscopic contact angle ($\theta_{\rm{s}}$) decreases as the distance
between pillars increases, therefore we expect $\theta_{\rm{s}}$ to approach
$\theta_{\rm{e}}$ as $\phi\to 0$.

####  3.1.3 Advancing contact angle

In this section, we develop a predictive equation for the advancing contact
angle of an interface in the Wenzel wetting state on a rough surface. Writing
the mechanical energy balance equation from part I

$\sum_{i<j}\sigma_{ij}\frac{\overrightarrow{A}_{ij}-\overleftarrow{A}_{ij}}{\sigma_{12}A_{\rm{CV}}}-\cos\theta_{\rm{meb}}-\overline{D}=0.$
(36)

Here, $\overrightarrow{A}_{ij}$ and $\overleftarrow{A}_{jj}$ are $ij$
interface areas entering and leaving the control volume respectively, as it
moves over the surface, $A_{\rm{CV}}$ is the projection of the area swept by
the control volume on a flat surface parallel to the average solid surface and
$\theta_{\rm{meb}}$ is the macroscopic angle of the interface. For the present
case - that is, Wenzel wetting state over a surface with square cross-
sectioned pillars arranged in a square array and the interface advancing in
surface periodicity direction - the various area terms can be expanded as

$\begin{split}\overrightarrow{A}_{2S}&=A_{\rm{CV}}+nA_{\rm{CV}}(4ah),\quad\overrightarrow{A}_{1S}=0,\quad\overrightarrow{A}_{12}=0,\\\
\overleftarrow{A}_{2S}&=0,\quad\overleftarrow{A}_{1S}=A_{\rm{CV}}+nA_{\rm{CV}}(4ah),\quad\overleftarrow{A}_{12}=0.\end{split}$
(37)

Here, $n$ is the number of pillars per unit area. Substituting the values from
equation (37) into equation (36) yields

$\cos\theta_{\rm{meb,a}}=(1+4\phi(h/a))\cos\theta_{\rm{e}}-\overline{D}.$ (38)

Here, $\theta_{\rm{meb,a}}$ is the advancing contact angle of the interface
calculated using the mechanical energy balance equation (36). From equations
(38), (32) and (34), the expression for advancing contact angle can hence be
written as

$\cos\theta_{\rm{meb,a}}=r\cos\theta_{\rm{e}}-\overline{D},$ (39)

where $r=(1+4\phi h/a$) is the Wenzel roughness ratio. Equation (39) only
depends upon the pillar aspect ratio ($h/a$), contact angle on the flat
surface ($\theta_{\rm{e}}$) and the surface area fraction of the pillars
($\phi$). The constants $A$, $B$ and $C$ can be predicted from the energy
minimization (example table 1), which makes equation (39) fully predictive in
nature.

###  3.2 Receding interface

In this section, we present the results for the receding motion of an
interface in the direction of surface periodicity and understand its dynamics,
and the dissipation in energy during the TPCL jumps. Referring back to figure
3, the advancing motion of fluid-1 is accompanied by the receding motion of
fluid-2 by the same amount. Therefore, we simulate the receding motion of an
interface with Young’s angle $\theta_{\rm{e}}$ by simulating the advancing
motion of the surrounding fluid with Young’s angle $180$°$-\theta_{\rm{e}}$.
Here we present the simulation results for the receding motion of fluid-1
($\theta_{\rm{e}}=72$°) by simulating the advancing of fluid-2 with
$\theta_{\rm{e}}=108$°. However, in the following discussion, we present the
results as if fluid-1 is receding with $\theta_{\rm{e}}=72$°.

####  3.2.1 Interface behaviour

Figure 18(a) shows the equilibrium morphologies of a pinned interface
projected on the symmetry plane as fluid-1 recedes ($\theta_{\rm{e}}=72$° and
$h/a=0.50$). Some of the equilibrium interface profiles before depinning are
plotted as dashed-dotted lines. The blue solid line represents the interface’s
first critical state, after which any further advancement will result in the
TPCL jump. The equilibrium interface profile after depinning (i.e. the second
critical state) is shown by the solid red line. In figure 18(b) we show the
equilibrium interface morphologies projected the symmetry plane for different
area fractions ($\phi=0.13,0.08,0.04$, $0.03$ and $\theta_{\rm{e}}=72$°,
$h/a=0.50$). The first and second critical states are shown by the solid and
dashed lines respectively.

Figure 18: Equilibrium morphologies during the receding motion of an interface
projected on the symmetry plane. The results are calculated for
$\theta_{\rm{e}}=72$° and pillar aspect ratio ($h/a$) of 0.5 (all the angles
are measured from fluid-1 side). (a) The interface can find an equilibrium
morphology while pinned at the pillar front and slides on the pillar top when
it recedes, finally reaching the first critical state at the rear face of the
pillar. Equilibrium morphologies before the interface reach the critical state
are shown by dashed-dotted lines. The corresponding macroscopic contact angle
values for the equilibrium interface profiles are also shown ($\phi=0.13$). At
$\theta_{\rm{s}}=41$°, the interface reaches the first critical state. The
projection of the equilibrium interface morphologies in the first and second
critical states are shown by solid blue and red lines respectively. (b)
Equilibrium interface morphologies projected on the symmetry plane capturing
the interfacial morphologies in the first (solid lines) and second critical
states (dashed lines), for area fractions, $\phi=0.13$ (red), $\phi=0.08$
(green), $\phi=0.04$ (blue), and $\phi=0.03$ (black) respectively.

A comparison between the equilibrium interface morphologies during the
advancing (figure 6(c)) and receding (figure 18(a)) motion of an interface
with $\theta_{\rm{e}}=72$° reveals that during the receding motion, the
interface can exist in an equilibrium state such that a portion of the TPCL
lies on the pillar tops. The TPCL slides on the pillar top before it finally
gets pinned at the rear face of the pillar where the interface exists in first
critical state before the TPCL jumps. However, the first critical interface
morphologies during advancing as well as the receding motion exist at the rear
face of the interface. Figure 19 shows the first critical state morphologies
of the interface for a range of area fractions, projected on the symmetry
plane as well as the domain wall ($y=W/2)$ for the receding motion of fluid-1
on a surface with $\theta_{\rm{e}}=72$° and pillar aspect ratio ($h/a$) of
0.35. Similar to the advancing interface case (figure 13), in the case of a
receding interface as well, we observe the interface curvature to increase
with the pillar area fraction.

Figure 19: Equilibrium interface morphologies just before depinning (first
critical state) as projected on the symmetry plane (shown by dashed-dotted
lines) and the domain wall ($y=W/2$, shown by solid lines) for a receding
interface on a surface with $\theta_{\rm{e}}=72$°, $h/a=0.35$ and area
fractions $\phi=0.01,0.03,0.04,0.08,0.13,0.20,0.30,0.40$ respectively. (b)
Zoomed view of the interface morphology near the TPCL.

####  3.2.2 Energy dissipation during receding motion of the interface

Figure 20 shows the variation in total non-dimensional energy dissipation
($\overline{D}$) and non-dimensional energy dissipation per pillar
($\overline{D}_{1}$) with the pillar area fraction ($\phi$) for a receding
interface on a surface with pillars of aspect ratio ($h/a$) 0.35 and 1.5 and
Young’s angle ($\theta_{\rm{e}}=72$°).

Figure 20: Variation in the total non-dimensional energy dissipation
($\overline{D}$, shown by triangles) and non-dimensional energy dissipation
per pillar ($\overline{D}_{1}$, shown by circles) with the pillar area
fraction ($\phi$) for an interface receding in the direction of surface
periodicity. Results have been simulated for Young’s angle $72$° and pillar
aspect ratio of (a) 0.35 and (b) 1.5 respectively. Fitting equations of dilute
(equation (32)) and non-dilute (equation (34)) form for $\overline{D}$ are
shown by dashed-dotted and solid black lines respectively. Fitting equations
for $\overline{D}_{1}$ are shown by dashed-dotted (dilute form, equation (30))
and solid (non-dilute, equation (33)) blue lines respectively.

The interface is moving in the direction of surface periodicity.
Qualitatively, the variation in $\overline{D}$ and $\overline{D}_{1}$ with
$\phi$ for a receding interface is similar to the variation in both the
dissipation ($\overline{D}_{1}$ and $\overline{D}$) for an advancing interface
(figure 17), which is expected as the receding motion of an interface can also
be identified as the advancing motion of the surrounding fluid. Therefore, we
use the fitting equations of the forms as given in equations (30), (32) for
$\overline{D}_{1}$ and $\overline{D}$ under the dilute surface assumption and
(33), (34) for $\overline{D}_{1}$ and $\overline{D}$ under non-dilute surface
assumption respectively. In table 2 we show the fitting coefficients and
$R^{2}$ value for the total non-dimensional dissipation according to equations
(32) and (34) for a receding interface with $\theta_{\rm{e}}=72$° and pillar
aspect ratios ($h/a$) 0.35 and 1.5 respectively. As per the advancing cases,
the two forms for $\overline{D}$ and $\overline{D}_{1}$, are able to
accurately fit the numerically generated results, with the non-dilute form
again producing slightly higher $R^{2}$ values.

Equation |  $h/a=0.35$ |  $h/a=1.5$
---|---|---
$\overline{D}=A\phi\ln\phi+C\phi$ (dilute surface) |  $A=-1.35$ $C=-0.10$ $R^{2}=0.975$ |  $A=-4.67$ $C=-6.85$ $R^{2}=0.859$
$\overline{D}=A\phi\ln\phi+B\phi^{2}+C\phi$ (non-dilute surface) |  $A=-2.09$ $B=3.54$ $C=-2.12$ $R^{2}=0.995$ |  $A=-7.46$ $B=38.78$ $C=-17.40$ $R^{2}=0.991$
Table 2: Fitting parameters for the total non-dimensional energy dissipation
($\overline{D}$) dependence on the pillar area fraction ($\phi$) for a
receding interface ($\theta_{\rm{e}}=72$°), as suggested in equations (32) and
(34).

From figures 17 and 20, we observe that the non-dimensional dissipation per
pillar ($\overline{D}_{1}$) decreases while the total non-dimensional
dissipation ($\overline{D}$) increases with the pillar area fraction ($\phi$)
for both the advancing and receding interfaces. Also, both of
$\overline{D}_{1}$ and $\overline{D}$ increase with the pillar aspect ratio
($h/a$) for advancing as well as receding interfaces. However, the magnitude
of non-dimensional dissipation per pillar and total non-dimensional
dissipation is different for the advancing and receding motion, even at the
same pillar area fraction. This is shown in figure 21 where we plot the
variation in $\overline{D}_{1}$ and $\overline{D}$ with $\phi$ for advancing
and receding motion of the interface on a surface with $\theta_{\rm{e}}=72$°
and pillar aspect ratios 0.35 and 1.5 respectively. We observe that the non-
dimensional dissipation per pillar and total non-dimensional dissipation is
higher for the advancing interface as compared to its receding motion on the
same surface. The energy dissipation is, therefore, asymmetric in nature which
may result in different rates at which the advancing contact angle increases
or the receding contact angle decreases with an increase in the pillar area
fraction.

Figure 21: Variation in the non-dimensional energy dissipation per pillar
($\overline{D}_{1}$, shown by circles) and total non-dimensional energy
dissipation ($\overline{D}$, shown by triangles) with the pillar area fraction
($\phi$) for an interface advancing (shown by filled symbols) and receding
(shown by empty symbols) in the direction of surface periodicity on a surface
with $\theta_{\rm{e}}=72$° and pillar aspect ratios ($h/a$) 0.35 (a) and 1.5
(b) respectively. The black lines are for equation (34) and blue lines for
equation (33) with solid and dashed representing the advancing and receding
motion of the interface.

####  3.2.3 Receding contact angle

Since the receding motion of fluid-1 is equivalent to the advancing motion of
fluid-2, the equation for the receding contact angle ($\theta_{\rm{meb,r}}$)
can be obtained from the equation for the advancing contact angle (equation
(39)) by replacing $\theta_{\rm{meb,a}}$ with (180°$-\theta_{\rm{meb,r}}$) and
$\theta_{\rm{e}}$ with (180°$-\theta_{\rm{e}}$) in equation (39), i.e.

$\begin{split}\cos\theta_{\rm{meb,r}}&=r\cos\theta_{\rm{e}}+\overline{D},\end{split}$
(40)

where $r=(1+4\phi h/a$) is the Wenzel roughness ratio. The total non-
dimensional energy dissipation ($\overline{D}$) is as per the form presented
in equations (32) and (34) with the coefficients as given in table 2.

###  3.3 Relationship between $\theta_{\rm{s}}$,
$\theta_{\rm{s}^{{}^{\prime}}}$ and $\theta_{\rm{meb}}$

Based on the results of our interface dynamics simulations, we propose three
main parameters for characterizing a typical interface advancement on a
structured surface in the direction of surface periodicity: the (first)
critical angle ($\theta_{\rm{s}}$), the macroscopic contact angle after the
TPCL depinning i.e. the second critical angle
($\theta_{\rm{s}^{{}^{\prime}}}$) and the dissipation in energy
($\overline{D}_{1}$ or $\overline{D}$). From figure 22 we observed that the
critical angle ($\theta_{\rm{s}}$) is independent of domain height provided
that the domain is high enough so as to ensure an approximately planar
interface near the top. However, $\theta_{\rm{s}^{{}^{\prime}}}$ depends upon
the domain height, which is shown in figure 22 where we plot the variation in
$\theta_{\rm{s}^{{}^{\prime}}}$ with $\phi$ for different domain heights ($H$)
expressed as a function of the domain width ($W$).

Figure 22: Variation in the macroscopic contact angle corresponding to the
second critical state ($\theta_{\rm{s^{{}^{\prime}}}}$) with the pillar area
fraction ($\phi$) for $\theta_{\rm{e}}=72$° and $h/a=0.5$. The (first)
critical contact angles ($\theta_{\rm{s}}$) are also shown as black circles.

We observe that $\theta_{\rm{s}^{{}^{\prime}}}$ increases as the domain height
is increased and approaches the (first) critical angle ($\theta_{\rm{s}}$,
shown as black circles) as the domain height becomes large. Therefore, we
expect that as the domain height approaches infinity, the interface should
advance with a constant macroscopic angle ($\theta_{\rm{s}}$). Also, since the
mechanical energy balance framework discussed in part I relates the
dissipation in energy during interface advancement to a single macroscopic
contact angle, we expect the advancing contact angle based on the balance of
mechanical energy (equation (39)) to approach $\theta_{\rm{s}}$ under the
assumption of separation of length scales i.e. large domain heights.

In figure 23 we plot the variation in $\theta_{\rm{meb,a}}$ with the pillar
area fraction ($\phi$) for different domain heights. For calculating the
advancing angle, we have used the non-dilute form of total non-dimensional
dissipation (equation (34)) in equation (39) as it represents the numerical
data most accurately. We observe that $\theta_{\rm{meb,a}}$ increases with the
domain height and approaches $\theta_{\rm{s}}$ (shown as black circles) at
higher domain heights.

Figure 23: Variation in the advancing contact angle ($\theta_{\rm{meb,a}}$)
calculated using equation (39) with total non-dimensional energy dissipation
calculated from energy minimization simulations under a non-dilute surface
assumption (equation (34)) with pillar area fraction ($\phi$) for different
domain heights on a surface with $\theta_{\rm{e}}=72$° and $h/a=0.5$. The
(first) critical contact angles ($\theta_{\rm{s}}$) are shown as black
circles.

The dependency of $\theta_{\rm{meb,a}}$ upon domain height originates from the
dependency of the energy dissipation upon the domain height. We observe that
the energy dissipation varies with the domain height and converges to a single
value as the domain height is increased (refer to Appendix F and figure 38).
From figures 22 and 23 we can therefore conclude that as $H\to\infty$,
$\theta_{\rm{meb,a}}\to\theta_{\rm{s}}$ and also
$\theta_{\rm{s}^{{}^{\prime}}}\to\theta_{\rm{s}}$. Since the receding motion
of a fluid is equivalent to the advancing motion of the surrounding fluid,
therefore, the above discussion also holds for the receding mode. Hence, via
observation, we have established a relationship between interface dynamics and
energy conservation/dissipation for this particular periodic surface topology.

###  3.4 Contact angle hysteresis

In figure 24(a) we plot the (first) critical advancing ($\theta_{\rm{s,a}}$)
and receding ($\theta_{\rm{s,r}}$) contact angles from the energy minimisation
simulations with pillar area fraction ($\phi$) for $\theta_{\rm{e}}=72$° and
pillar aspect ratios ($h/a$) of 0.35 and 1.5 respectively, along with the
mechanical energy balance derived values ($\theta_{\rm{meb,a}}$,
$\theta_{\rm{meb,r}}$) determined via equations (39) and (40). The mechanical
energy balance values use the non-dilute form of total non-dimensional energy
dissipation for greater accuracy. Under these conditions, the critical and
mechanical energy balance derived angles are very similar across all tests
(refer to figure 23). We observe that the advancing contact angle increases
with the pillar area fraction for both pillar aspect ratios, with higher
advancing contact angles being achieved at higher aspect ratios for any value
of $\phi$. Similarly, the receding contact angle decreases with increases in
$\phi$, with lower value angles being found at higher pillar aspect ratios. We
also observe that the increase in $\theta_{\rm{meb,a}}$ with $\phi$ is greater
than the decrease of $\theta_{\rm{meb,r}}$ with $\phi$. This asymmetry in the
contact angles is due to the asymmetry in the energy dissipation between the
advancing and receding modes.

Figure 24: Variation in the advancing and receding contact angles obtained
using energy minimization simulations ($\theta_{\rm{s,a}}$,
$\theta_{\rm{s,r}}$) and mechanical energy balance equations (39) and (40)
($\theta_{\rm{meb,a}}$, $\theta_{\rm{meb,r}}$) with the pillar area fraction
($\phi$) for pillar aspect ratios ($h/a$) 0.35 and 1.5 respectively and
Young’s angle, $\theta_{\rm{e}}=72$° (a) and 108° (b) respectively. (c) Shows
the variation in contact angle hysteresis in degrees, obtained using energy
minimization simulations (i.e.
$\Delta\theta_{\rm{s}}=\theta_{\rm{s,a}}-\theta_{\rm{s,r}}$) and mechanical
energy balance equations (39) and (40) (i.e.,
$\Delta\theta_{\rm{meb}}=\theta_{\rm{meb,a}}-\theta_{\rm{meb,r}}$) with the
pillar area fraction ($\phi$) for pillar aspect ratios $h/a=0.35$ and 1.5
respectively and $\theta_{\rm{e}}=72$°.

In figure 24(b) we plot $\theta_{\rm{s,a}}$, $\theta_{\rm{s,r}}$,
$\theta_{\rm{meb,a}}$ and $\theta_{\rm{meb,r}}$ with $\phi$, for
$\theta_{\rm{e}}=108$° and pillar aspect ratios ($h/a$) of 0.35 and 1.5
respectively. Qualitatively, the contact angle variation with $\phi$ for
$\theta_{\rm{e}}=72$° and $\theta_{\rm{e}}=108$° are similar, but for
$\theta_{\rm{e}}=108$° we observe that the decrease in $\theta_{\rm{meb,r}}$
with $\phi$ is greater than the increase in $\theta_{\rm{meb,a}}$ with $\phi$.
Since the receding motion of fluid-1 is equivalent to the advancing motion of
fluid-2, therefore, the receding angle variations for $\theta_{\rm{e}}=72$°
are equivalent to the variations in the advancing angle for
$\theta_{\rm{e}}=108$°. Similarly, the advancing angle variation for
$\theta_{e}=72$° is equivalent to the receding angle variations for
$\theta_{\rm{e}}=108$°.

In figure 24(c) we plot the variation in contact angle hysteresis
($\Delta\theta_{\rm{s}}=\theta_{\rm{s,a}}-\theta_{\rm{s,r}}$ and
$\Delta\theta_{\rm{meb}}=\theta_{\rm{meb,a}}-\theta_{\rm{meb,r}}$) with the
pillar area fraction ($\phi$) for pillar aspect ratios ($h/a$) of 0.35 and 1.5
respectively and $\theta_{\rm{e}}=72$°. The contact angle hysteresis based on
the mechanical energy balance is calculated using equations (39) and (40) with
the non-dilute form of total non-dimensional dissipation ($\overline{D}$). As
expected, the contact angle hysteresis increases with both the pillar area
fraction ($\phi$) and the pillar aspect ratio ($h/a$). However, as the area
fraction is reduced, the hysteresis also decreases and eventually approaches
zero as $\phi\to 0$. This is due to zero total non-dimensional energy
dissipation, as the surface becomes very dilute (equation (35)). In reality,
there may be non-zero hysteresis even when $\phi\to 0$ due to the presence of
roughness at smaller length scales, or other causes of inherent hysteresis.

###  3.5 Comparison with experiments

In this section we compare our numerical model with the experimental data of
Forsberg et al. (2010). Forsberg et al. measured advancing and receding
contact angles on surfaces decorated with microscopic pillars of square cross-
section with a fixed width of 20 $\mu\rm{m}$ and heights 7 $\mu\rm{m}$ and 30
$\mu\rm{m}$ respectively. The advancing and receding contact angles measured
on a flat surface were $72$° and $59$° respectively. Before we proceed
further, it is important to address the presence of inherent hysteresis on
surfaces. In developing our model (part I) we assumed the surface to be
perfectly smooth at length scales below $h_{\rm{rough}}$, chemically
homogeneous and with no irreversibility in the creation or destruction of
surfaces on a molecular scale. The presence of any one of these effects may
result in a non-zero inherent hysteresis in the system. For example, all
surfaces contain roughness at different length scales irrespective of the
manufacturing methods used (Quéré, 2008), which can result in a non-zero
inherent hysteresis.

Figure 25: Total non-dimensional dissipation ($\overline{D}$) for receding
interface with $\theta_{\rm{R}}=59$° and pillar aspect ratios 0.35 and 1.5
respectively. Inset shows that at higher area fractions, we may observe a
permanently pinned TPCL.

In order to accommodate the inherent hysteresis on the surface, instead of
$\theta_{\rm{e}}$ we use $\theta_{\rm{A}}$ and $\theta_{\rm{R}}$ for advancing
and receding motion of the interface in our analysis (He et al., 2004), where
$\theta_{\rm{A}}$ and $\theta_{\rm{R}}$ are the advancing and receding contact
angles measured on the flat surface. Substituting $\theta_{\rm{e}}$ with
$\theta_{\rm{A}}$ and $\theta_{\rm{R}}$ in equations (39) and (40)
respectively, gives

$\cos\theta_{\rm{meb,a}}=r\cos\theta_{\rm{A}}-\overline{D},$ (41)
$\cos\theta_{\rm{meb,r}}=r\cos\theta_{\rm{R}}+\overline{D}.$ (42)

The total non-dimensional energy dissipation for an advancing interface with
$\theta_{\rm{A}}=72$° and pillar aspect ratios $0.35$ and $1.5$ is shown in
figure 17 and the values of coefficients in the dissipation equation (34) are
given in table 1. In figure 25 we show the variation of total non-dimensional
dissipation ($\overline{D}$) with $\phi$ for the receding interface with
$\theta_{\rm{R}}=59$° and pillar aspect ratios $0.35$ and $1.5$ respectively.
We use equations (41) and (42) and the dissipation parameters as tabulated in
table 3 to calculate the advancing and receding contact angles.

System |  $h/a=0.35$ |  $h/a=1.5$
---|---|---
Advancing interface $\theta_{\rm{A}}=72$° |  $A=-1.80$ $B=1.09$ $C=0.77$ |  $A=-11.26$ $B=53.16$ $C=-20.95$
Receding interface $\theta_{\rm{R}}=59$° |  $A=-1.83$ $B=4.72$ $C=-2.73$ |  $A=-7.04$ $B=2.79$ $C=-16.33$
Table 3: Fitting parameters for the total non-dimensional energy dissipation
($\overline{D}$) dependence on the pillar area fraction ($\phi$) for advancing
and receding interface on a surface with $\theta_{\rm{A}}=72$°,
$\theta_{\rm{R}}=59$° and pillar aspect ratio of 0.35 and 1.5 respectively.
The form of the fitting equation is as suggested in equation (34). Figure 26:
Variation in the advancing ($\theta_{\rm{a}}$) and receding
($\theta_{\rm{r}}$) contact angles (in degrees) with the pillar area fraction
($\phi$) for Wenzel wetting state. The experimental results of Forsberg et al.
(2010) ($\theta_{\rm{exp,a}}$, $\theta_{\rm{exp,r}}$), numerical simulations
($\theta_{\rm{s,a}}$, $\theta_{\rm{s,r}}$) and the mechanical energy balance
equation (41), (42) ($\theta_{\rm{meb,a}}$, $\theta_{\rm{meb,r}}$) are shown
for a surface with $\theta_{\rm{A}}=72$°, $\theta_{\rm{R}}=59$° and pillar
aspect ratio of 0.35 and 1.5 respectively. Wenzel’s equation (2) for the
advancing and receding interface is also plotted for comparison. The advancing
contact angle based on the numerical simulations ($\theta_{\rm{s45,a}}$) and
the balance of mechanical energy using equations (41)
($\theta_{\rm{meb45,a}}$) for an interface moving at an angle of 45° to the
surface periodicity direction ($\theta_{\rm{A}}=72$°, $h/a=0.35$) is also
shown.

Figure 26 shows experimental data from Forsberg et al. (2010) along with our
simulation results from equations (41) and (42). The (first) critical contact
angles ($\theta_{\rm{s}}$) obtained from the simulations for both the
advancing and receding interface are also shown. Very good agreement between
the simulation results and experimental data is observed at low pillar area
fractions. At higher area fractions, the proposed equations tend to
overpredict the advancing contact angles slightly. Receding contact angle
values are in good agreement with our model for both the aspect ratios
($h/a=0.35$ and 1.5) for the entire range of experimental data, which is up to
$\phi=0.21$ for $h/a=0.35$ and $\phi=0.04$ for $h/a=1.5$. It should be noted
that at higher area fractions the receding TPCL can get permanently pinned at
the start on the first pillar, causing the receding contact angle to approach
$0$° and a film of liquid to be left within the pillars (see figure 25 inset).
The onset of an equilibrium morphology with a permanently pinned TPCL depends
upon the pillar geometry and the receding contact angle on a flat surface
($\theta_{\rm{R}}$). For example on a surface with $h/a=1.5$ and
$\theta_{\rm{R}}=59$° we observe that for $\phi\geq 0.13$ the TPCL gets
permanently pinned. The morphological transitions in such cases have already
been discussed in the context of an advancing interface (see figure 16), and
even though the present numerical method is capable of dealing with these
cases study of this transition is left to future work. Also, at a certain low
Young’s angle the TPCL can form a closed loop around the pillars (Semprebon et
al., 2012) such that a portion of the pillar is submerged under the liquid or
the liquid may start exhibiting hemiwicking (Ishino et al., 2007). The study
of such equilibrium interfacial morphologies is also left to future work.

To understand the small differences between the simulated and experimental
results shown in figure 26, particularly for the advancing state, it is
important to note that during the spreading of a droplet, the interface can
approach a pillar from different directions. Figure 27 (inset) shows a TPCL
approaching pillars at an angle $\psi$ relative to the periodicity direction
of the structured surface.

Figure 27: Representation of the TPCL on a surface with a structured array of
pillars. The TPCL can approach pillars from different directions (a typical
scenario is shown in the inset where the TPCL is approaching pillars at an
angle of $\psi$°), resulting in different equilibrium morphologies of the
interface. Inset shows interface morphologies during the critical state for
different approach ($\psi$) directions. Interface morphologies are simulated
for a surface with $\theta_{\rm{e}}=72$° and $h/a=0.35$. The interface shapes
in the inset only capture the effect of the pillar shape when the contact line
approaches from a certain direction. In the actual scenario, the pillar shape
as well as the inter-pillar distance would change depending on the direction
of the approaching interface and we may end up getting very different
interface morphologies.

In order to show the impact of the interface movement direction with respect
to the surface periodicity, in figure 26 we also plot the critical contact
angle ($\theta_{\rm{s45,a}}$) and the mechanical energy balance derived
contact angle ($\theta_{\rm{meb45,a}}$) for an interface advancing at an angle
($\psi$) of 45° ($\theta_{\rm{e}}=72$° and $h/a=0.35$) that we have computed
numerically using a slightly different computational domain (see Appendix G).
We observe that $\theta_{\rm{s45,a}}$ values for $\psi=45$° are less than the
equivalent $\theta_{\rm{s,a}}$ values calculated for an interface advancing in
the direction of the surface periodicity. Given that the actual advancing or
receding contact angle of a droplet spreading on a structured surface would be
an average value for the interface advancing in different directions relative
to the surface periodicity (Dorrer & Rühe (2008) made a similar suggestion) it
is encouraging to see that the experimentally measured values for
$\theta_{\rm{exp,a}}$ are bounded by the equivalent $\theta_{\rm{meb,a}}$ and
$\theta_{\rm{meb45,a}}$, results. In fact, analysis of the Forsberg et al.
(2010) study suggests that the advancing and receding contact angles were
measured in such a way that they capture the TPCL movement in the direction of
surface periodicity. Therefore it is encouraging that the
$\theta_{\rm{exp,a}}$ scatter leans towards the $\theta_{\rm{s,a}}$ values
rather than the $\theta_{\rm{s45,a}}$ values. We have shown the results for
$\psi=45$° since the surface is periodic along that direction and we can use
periodic boundary conditions in the simulations (as detailed in Appendix G).
However, for other directions of interface motion we cannot use periodic
boundary conditions and a much wider computational domain is required to
simulate the interface dynamics. These simulations, while feasible, are left
to future study.

##  4 Conclusions

We have presented a novel numerical method to simulate the microscale
interface dynamics on a structured surface. The present model can predict
contact angle hysteresis during homogeneous wetting (Wenzel state) of micro-
structured surfaces from just the knowledge of surface topology. By using
superquadrics to model the pillar geometry, we are able to simulate the
interface dynamics over pillars with rounded edges, closely. Using a surface
minimization approach to approximate fluid-driven interface movement, we
obtain the energy dissipation during the advancing and receding motion of the
interface.

With our simulations that utilize a constant Young’s angle, at low pillars
area fractions we observe a logarithmic dependence of total non-dimensional
energy dissipation ($\overline{D}$) and non-dimensional energy dissipation per
pillar ($\overline{D}_{1}$) on the pillar area fraction ($\phi$), broadly
consistent with the Joanny & de Gennes (1984) analysis of pinning on sparsely
spaced (dilute) strong defects. At higher area fractions a correction term
($B\phi$) was added to capture the predicted dissipation. Both
$\overline{D}_{1}$ and $\overline{D}$ depend upon the pillar aspect ratio
($h/a$) as well, with higher dissipation observed at higher aspect ratios. We
simulated the receding motion of an interface with Young’s angle
($\theta_{\rm{e}}$) by simulating the advancing motion of surrounding fluid
with Young’s angle (180°- $\theta_{\rm{e}}$). We observed a logarithmic
relationship between $\overline{D}_{1}$ & $\phi$ and $\overline{D}$ & $\phi$
for the receding interface as well. The dissipation for the receding interface
also increases with the pillar aspect ratio. At certain higher aspect ratios
and/or higher area fractions we observed permanent pinning of the TPCL
indicating the maximum range of validity of our Wenzel-based results. We used
the roughness scale mechanical energy balance from part I to obtain the
equation for advancing and receding contact angles on surfaces with no
inherent hysteresis (equations (39), (40)) and on surfaces with a finite
inherent hysteresis (equations (41), (42)), based on the results of energy
dissipation calculations.

Finally, we compared our model with the experimental data of Forsberg et al.
(2010) and obtained very good agreement between simulation and experimental
results across all relevant area fractions, but particularly at low area
fractions. At higher area fractions we find the model slightly overpredicts
the advancing contact angle. However, we also observed that the maximum stable
contact angle ($\theta_{\rm{s}}$) for the interface depends upon the direction
from which the interface approaches a pillar. Isolated simulations using an
advance direction of ($\psi=45$°) to the surface periodicity direction show
that neither the $\psi=0$° or 45° is the best representation of the interface
advancement, but that the experimental results are in fact ‘bracketed’ by the
simulation results at these two angles. We conclude that, orientation effects
should be more fully considered in future work. In this work, we have
considered only a square pillar geometry using a very limited number of
Young’s angles and pillar aspect ratios. It would be interesting to see the
effect of different pillar geometries and Young’s angle on contact angle
hysteresis. Again, this is left to future work. Lastly, the micro-scale
interface dynamics and the CAH calculated on pillared surfaces are observed to
depend upon the pillar aspect ratio rather than their actual size. Therefore,
the present model for simulating interface advancement and the dissipation
equations presented in §3.1.2 are valid for any pillar sizes for which the
separation of length scales can be observed in the system. Even though we
developed equations (30) to (34) for the roughness in the form of micron sized
pillars, the same equations can be used on surfaces with nanometric roughness
as well, provided that the defects are strong enough to pin the interface, and
that the various velocity and length scale constraints detailed in part I are
satisfied.

## Acknowledgements

One of the authors (P.K.) acknowledges the financial support from the
University of Melbourne in form of the Melbourne research Scholarships
program. P.K. also acknowledges the support from Melbourne India Postgraduate
Program (MIPP).

## Declaration of interests

The authors report no conflict of interest.

## Author ORCIDs

Pawan Kumar, https://orcid.org/0000-0002-8654-8514
D.J.E. Harvie, https://orcid.org/0000-0002-8501-1344

##  Appendix A Effect of step size on the interface dynamics

In §2.2.4 we discussed the method of incremental advance for simulating the
interface dynamics on a rough surface. There we commented that a small $\Delta
x$ is desired for accurately capturing the interface movement. In figure 28 we
show the variation in macroscopic contact angle ($\theta_{\rm{m}}$) with the
interface top position ($x_{\rm{top}}$) as a function of $\Delta x$.
Specifically we have used seven different step sizes ($\Delta
x=0.10,0.25,0.50,1.0,2.0,4.0$ and 8.0 respectively) for simulating the
interface advancement using the method of incremental advance on a surface
with $\phi=0.08$, $h/a=0.5$ and $\theta_{\rm{e}}=72$° within a domain which is
$3.5$ units wide and $35$ units high. Even though we observe the same critical
angle ($\theta_{\rm{s}}$) for $\Delta x=1.0,2.0,4.0$ and 8.0, this, however,
may not be the most accurate prediction of $\theta_{\rm{s}}$. Since the
accuracy with which $\theta_{\rm{s}}$ is calculated depends upon the magnitude
of the step size, we observe a different $\theta_{\rm{s}}$ when a smaller step
size is used. For example, we show the $\theta_{\rm{m}}$ vs. $x_{\rm{top}}$
variation capturing only the first and second critical states of the interface
when step sizes of $\Delta x=0.50,0.25$ and 0.10 are used (shown as inset). We
can see that the $\theta_{\rm{s}}$ value depicting the critical state of the
interface depends upon the magnitude of $\Delta x$ with a better prediction of
the critical state for a smaller $\Delta x$. Another observation from the
figure is that for $\Delta x=4.0$, $\theta_{\rm{s}^{{}^{\prime}}}$ is greater
than $\theta_{\rm{s}}$, which is due to an inaccurate prediction of
$\theta_{\rm{s}}$ when a bigger step size ($\Delta x=4.0$) is used.

Another interesting result that we may observe when a bigger step size is used
for interface advancement is that the interface may not be able to find an
equilibrium position on the next pillar in the direction of interface
advancement and may skip a few pillars before the TPCL can finally pin again
at a suitable location. We observed such behaviour for a step size of $\Delta
x=8.0$, where the TPCL after depinning from the first pillar skipped two
pillars in the direction of advancement before pinning again at the third
pillar. Equilibrium interface morphologies depicting this behaviour are shown
in figure 29.

Figure 28: Variation in the macroscopic contact angle ($\theta_{\rm{m}}$) with
the position of the interface ($x_{\rm{top}}$) for an advancing interface
simulated by incremental advance method with a step size of $\Delta
x=1.0,2.0,4.0$ and 8.0 respectively. The inset shows a zoomed view of the
variation in $\theta_{\rm{m}}$ with $x_{\rm{top}}$ around a first critical
state which also shows the variation with three additional step sizes of 0.50,
0.25 and 0.10 respectively, capturing only the first and second critical
states. The simulations are performed with $\theta_{\rm{e}}=72$°, $h/a=0.35$,
$\phi=0.08$, $W=3.5$ and $H=10W$. Figure 29: (a) Equilibrium interface
morphologies capturing the TPCL jump during the interface advancement for two
different step sizes, $\Delta x=4.0$ and $\Delta x=8.0$. The equilibrium
interface morphology representing the second critical state for $\Delta x=8.0$
is shown red in colour. For $\Delta x=8.0$, the TPCL skips two pillars in the
advancement direction before getting pinned on the third pillar. (b)
Equilibrium interface morphologies depicting TPCL jump for $\Delta x=4.0$ and
8.0 respectively, showing the full domain height. The simulations are
performed with $\theta_{\rm{e}}=72$°, $h/a=0.35$, $\phi=0.08$, $W=3.5$ and
$H=10W$.

##  Appendix B Equilibrium interface morphologies during a typical advancement

Equilibrium interface morphologies depicting a typical interface advancement
are shown in figure 30. The simulation run starts with an equilibrium
morphology with $x_{\rm{top}}$ chosen such that
$\theta_{\rm{m}}=\theta_{\rm{e}}$. As the interface is advanced, the TPCL
remains pinned on the pillar until the first critical state is reached
($\theta_{\rm{m}}=\theta_{\rm{s}}$). Upon further advancement of the
interface, the TPCL depins from the first pillar and the interface only finds
an equilibrium morphology (second critical state) when the TPCL gets pinned on
the next pillar in the direction of advancement
($\theta_{\rm{m}}=\theta_{\rm{s}^{{}^{\prime}}}$). The TPCL remains pinned on
this pillar until the first critical state is reached again upon further
interface advancement (i.e. $\theta_{\rm{m}}$ increases from
$\theta_{\rm{s}^{{}^{\prime}}}$ to $\theta_{\rm{s}}$) at which it executes a
jumping event. In figure 30 we have shown three TPCL jumping events depicting
its typical stick-slip motion.

Figure 30: Equilibrium interface morphologies depicting a typical advancing
motion. The TPCL moves in a stick-slip manner executing jumps, three such
jumping events are shown in the figure. The simulations are performed with
$\theta_{\rm{e}}=72$°, $h/a=0.35$, $\phi=0.08$, $W=3.5$ and $H=10W$.

##  Appendix C The critical-state method for simulating interfacial dynamics

In this section, we explain the critical-state method for simulating the
interface dynamics on a rough surface. This method involves two different
steps for obtaining the first and second critical states respectively. For the
first critical state, we start with an initial equilibrium morphology obtained
by minimizing the total energy within the domain consisting of a flat surface
inclined at Young’s angle relative to the domain base. To carry out energy
minimization, the boundary conditions used are
$\theta_{\rm{m}}=\theta_{\rm{e}}$ at the domain top, equation (20) along the
TPCL and the orthogonality between the interface and the domain walls
(equation (19)). If $d$ is the inter-pillar distance and the center of the
pillar is situated at $y=0$, the domain walls are situated at $y=-d/2$ and
$y=d/2$ respectively, on planes parallel to the $x-z$ plane (and at $y=0$ and
$y=d/2$ if a half domain is used, see figure 3 for domain geometry and
coordinate system). Also, initially, the interface is located near the pillar
such that it is intersecting it. The intersecting facets of the pillar and the
interface are identified and the intersecting facets, edges and vertices (of
the interface) are constrained to follow the shape of the pillar, i.e. they
are constrained to follow equation (18). Since the portion of the interface
constrained to follow the pillar’s profile represents the fluid-1/solid
interface instead of the fluid-1/fluid-2 interface, its surface energy is
changed from $\sigma_{12}$ to $-\sigma_{12}\cos\theta_{\rm{e}}$. Now the
energy minimization is carried out in a number of steps with gradual mesh
refinements and an equilibrium interface morphology is obtained. This is the
initial equilibrium shape with $\theta_{\rm{m}}=\theta_{\rm{e}}$. To simulate
the advancing motion of the interface, $\theta_{\rm{m}}$ is increased by a
small amount ($\Delta\theta$), which is the step size of interface advancement
and the energy minimization with gradual refinements is carried out again to
obtain the equilibrium interface morphology with
$\theta_{\rm{m}}=\theta_{\rm{m}}+\Delta\theta$, with the rest of the boundary
conditions unchanged. This process is repeated to simulate the interface
advancement. The interface advancement is stopped when the $\theta_{\rm{m}}$
reaches a value such that any further advancement results in the depinning of
the TPCL, i.e. the interface cannot exist in equilibrium with the TPCL pinned
on the pillar. This equilibrium interface morphology is the first critical
morphology and the macroscopic angle is the (first) critical angle
($\theta_{\rm{s}}$) (see Forsberg et al. (2010), Semprebon et al. (2012) for a
discussion on the advancing contact angle being the same as the maximum
macroscopic contact angle for which the interface can be in mechanical
equilibrium while pinned on a pillar). The accurate prediction of the first
critical interface morphology depends upon the step size and a smaller step
size gives a better approximation. Therefore, once a first critical state is
obtained by using a certain $\Delta\theta$, the step size is reduced to
$\Delta\theta/2$ and the first critical interface morphology is obtained
again. To this end, the first critical interface morphology obtained using the
step size of $\Delta\theta$ serves as the starting point and the process of
successive energy minimization and interface advancement is repeated to obtain
a better estimate of the first critical interface morphology. This process is
repeated again until we reach the minimum step size, $\Delta\theta_{\rm{min}}$
(for example, 0.1° which is used in this study). The discussion here covers
the first step of the method, which is obtaining the first critical interface
morphology. Figure 31(a) shows the algorithm for this step of the method.

Once the first critical interfacial morphology is obtained to the desired
level of accuracy, we proceed to obtain the second critical state of the
interface. Based on the assumption of a slow-moving interface, the second
critical state is the equilibrium interface morphology calculated such that
the interface top remains at the same location, i.e. $x_{\rm{top}}=x_{\rm{s}}$
($x_{\rm{s}}$ is obtained from the first step) while the TPCL pins on the next
pillar along the advancing direction. To obtain the second critical state we
start with an initial equilibrium interface morphology which is obtained by
minimizing the total energy within the domain containing a flat interface
inclined at Young’s angle relative to the domain bottom and positioned close
to the second pillar in the advancement direction such that the interface is
touching the pillar. The boundary conditions used here are
$x_{\rm{top}}=x_{\rm{initial}}$ at the top, equation (20) along the TPCL and
equation (19) along the domain walls. Here, $x_{\rm{initial}}$ is the initial
value for positioning the interface such that
$x_{\rm{initial}}=x_{0}-H/\tan\theta_{\rm{e}}$, this ensures that we start
with macroscopic contact angle being same as Young’s angle. The intersecting
facets of the second pillar and the interface are identified and the facets,
edges and vertices (of the interface) are constrained to follow the profile of
the intersecting pillar. Also, the surface energy of the intersecting facets
of the interface is changed to $-\sigma_{\rm{12}}\cos\theta_{\rm{e}}$ (from
$\sigma_{12}$) followed by energy minimization with gradual mesh refinements.
The equilibrium interface morphology thus obtained is the initial morphology
which is used to arrive at the second critical morphology by gradual
advancements and energy minimization (with gradual mesh refinements). The
interface is advanced by moving the top, i.e.
$x_{\rm{top}}=x_{\rm{top}}+\Delta x$, where $\Delta x$ is the step size. The
interface advancement is carried out until $x_{\rm{top}}+\Delta x<x_{\rm{s}}$,
and when this happens, $x_{\rm{top}}$ is set as equal to $x_{\rm{s}}$ followed
by the energy minimization with gradual mesh refinements. The equilibrium
interface morphology thus obtained (at $x_{\rm{top}}=x_{\rm{s}}$) is the
second critical morphology. This concludes the second and final step of the
process. A flow chart depicting the algorithm used for obtaining the second
critical state is shown in figure 31(b).

(a) Step1 - Obtaining the first critical state of the interface.

(b) Step2 - Obtaining the second critical state of the interface.

Figure 31: Flowchart of the algorithm used for simulating interfacial dynamics
using the critical-state method.

In figure 32 we show a typical variation in the total non-dimensional energy
within the simulation domain ($\overline{E}$)333When the interface top is
constrained by an angle ($\theta_{\rm{m}}$), additional energy equal to
$x_{\rm{top}}W\cos\theta_{\rm{m}}$ is added to the total energy by SE, which
needs to be subtracted from the total SE energy to get the true interfacial
energy within the simulation domain. and the macroscopic contact angle
($\theta_{\rm{m}}$) with the interface position ($x_{\rm{top}}$) during the
interface advancement as captured by the critical-state method. Simulation
results obtained with two different step sizes ($\Delta\theta$) i.e.,
$\Delta\theta=2.0$° (filled symbols) $\Delta\theta=1.0$° (empty symbols) are
shown in the figure. We find that irrespective of the step size used, a very
similar variation in the total non-dimensional energy and the macroscopic
angle is observed. However, the critical states identified by the method can
be different for different step sizes used (i.e., $\Delta\theta$, or more
specifically $\Delta\theta_{\rm{min}}$ used in the simulation). For example,
in figure 32 we have used two different $\Delta\theta_{\rm{min}}$, 0.1° for
the simulation run with $\Delta\theta=2.0$° and 1.0° for the run with
$\Delta\theta=1.0$° respectively. We can observe two different first critical
states (two different $\theta_{\rm{s}}$) identified by the method for two
different $\Delta\theta_{\rm{min}}$ values, with the better approximation of
the critical-state being achieved at the smaller $\Delta\theta_{\rm{min}}$.
Since the second critical state depends upon the interface position during the
first critical state ($x_{\rm{s}}$), a smaller $\Delta\theta_{\rm{min}}$
ensures a better approximation of the second critical state as well.

Figure 32: Variation in the total non-dimensional energy within the simulation
domain ($\overline{E}$, shown as blue circles) and macroscopic contact angle
($\theta_{\rm{m}}$, shown as black squares) with the position of the interface
($x_{\rm{top}}$) for an advancing interface simulated by critical-state method
with step sizes of $\Delta\theta=2.0$° (shown by filled symbols) and
$\Delta\theta=1.0$° (shown by empty symbols). For $\Delta\theta=2.0$° we use
$\Delta\theta_{\rm{min}}=0.1$°. The simulations are performed with
$\theta_{\rm{e}}=72$°, $h/a=0.35$, $\phi=0.08$, $W=3.5$ and $H=10W$.

##  Appendix D Comparison between the incremental advance and critical-state
method

In this work, we have demonstrated two methods, namely the incremental advance
method and the critical-state method for simulating the advancing motion of an
interface. Here we demonstrate that both methods capture the same interfacial
dynamics. In figure 33 we show the variation in total non-dimensional energy
within the simulation domain ($\overline{E}$) and the macroscopic contact
angle ($\theta_{\rm{m}}$) with the interface position ($x_{\rm{top}}$) for an
advancing interface simulated by incremental advance (filled symbols) and
critical-state (empty symbols) methods respectively. The step sizes used for
simulating the interface advancement are $\Delta x=2.0$ for the incremental
advance method and $\Delta\theta=2.0$°, $\Delta\theta_{\rm{min}}=0.1$° for the
critical-state method. We observe a very similar variation in $\overline{E}$
and $\theta_{\rm{m}}$ with $x_{\rm{top}}$ predicted via both the methods. The
only difference is in the precise identification of the critical states, which
is due to the finite step size used in the simulations. For example, when a
smaller $\Delta x$ is used in the incremental advance method, the first
critical state identified by the incremental advance method approaches that
identified by the critical-state method. In figure 33 we show the variation in
$\overline{E}$ and $\theta_{\rm{m}}$ with $x_{\rm{top}}$ capturing the first
and second critical states when a step size of $\Delta x=0.10$ is used in the
incremental advance method as red circles and squares respectively. We observe
that the first and second critical states as captured by the incremental
advance method approach to the corresponding states captured by the critical-
state method when a smaller $\Delta x$ (0.1) is used. We therefore expect that
this small difference between the two methods should vanish when a very small
step size, i.e. $\Delta x\approx 0$ and $\Delta\theta_{\rm{min}}\approx 0$ is
used. The critical-state method is preferred over the incremental advance
method since we are only interested in the first and second critical states
for calculating the dissipation in energy during interface advancement and
hence this method is more efficient for a given accuracy level. However, there
is a limitation to this method in that that it can only be used on structured
surfaces when the interface is advancing in the direction of surface
periodicity and all the pillars are identical.

Figure 33: Variation in the total non-dimensional energy within the simulation
domain ($\overline{E}$, shown as blue circles) and macroscopic contact angle
($\theta_{\rm{m}}$, shown as black squares) with the interface position
($x_{\rm{top}}$) for an advancing interface simulated by incremental advance
method with a step size ($\Delta x$) of 2.0 (shown by filled symbols) and by
critical-state method with a step size ($\Delta\theta$) 2.0 (shown by empty
symbols). $\overline{E}$ and $\theta_{\rm{m}}$ values capturing the first and
second critical states for $\Delta x=0.1$ are also shown in red colour for the
incremental advance method. The simulations are performed with
$\theta_{\rm{e}}=72$°, $h/a=0.35$, $\phi=0.08$, $W=3.5$ and $H=10W$.

In the preceding discussion, we demonstrated that the incremental advance and
critical-state methods generate the same critical interface topologies as the
step sizes used in both methods approach zero. Therefore, we expect that the
energy dissipation calculated by the two methods should approach the same
value for small step sizes. In figure 34 we plot the variation in the non-
dimensional energy dissipation per pillar ($\overline{D}_{1}$) calculated by
the method of incremental advance with the step size ($\Delta x$). The
$\overline{D}_{1}$ value calculated by the critical-state method for a
$\Delta\theta_{\rm{min}}$ of 0.1° is also plotted (in red colour). We observe
that as the step size ($\Delta x$) is reduced, the non-dimensional dissipation
per pillar ($\overline{D}_{1}$) increases and approaches the
$\overline{D}_{1}$ calculated by the method of critical-state. Hence we
conclude that for an infinitesimal step size ($\Delta\theta\to 0$ and $\Delta
x\to 0$), the dissipation values calculated by either of the two methods will
be the same for an interface advancing on a structured surface in the
direction of surface periodicity.

Figure 34: Variation in the non-dimensional energy dissipation per pillar
($\overline{D}_{1}$) with the step size ($\Delta x$) calculated using the
incremental advance method for simulating interface dynamics.
$\overline{D}_{1}$ calculated using the critical-state method for a
$\Delta\theta_{\rm{min}}=0.1$° is a also shown as a red circle. The
simulations are performed with $\theta_{\rm{e}}=72$°, $h/a=0.35$, $\phi=0.08$,
$W=3.5$ and $H=10W$.

##  Appendix E Effect of mesh resolution

Surface Evolver uses the finite element discretization method with each
surface divided into a large number of triangular facets. The larger the
number of facets, the better the representation of the curved surface, but at
the same time the time required for convergence increases. We say the solution
has converged when there is minimal change in the total energy ($E$) of the
system with further iterations. Specifically, if $E^{n}_{g}$ represents the
total energy during the $n{\rm{th}}$ iteration then our convergence criteria
can be represented by

$|{E^{{n}}_{g}-E^{n+1}_{g}}|<\epsilon,$ (43)

where $\epsilon\ll 1$. We have used a value of $\epsilon=10^{-7}$. Since the