## Abstract

We develop a two-dimensional model for the transient diffusion of gas from the cavities in ridge-type structured surfaces to a quiescent liquid suspended above them in the Cassie state to predict the location of the liquid vapor-interface (meniscus) as a function of time. The transient diffusion equation is numerically solved by a Chebyshev collocation (spectral) method coupled to the Young–Laplace equation and the ideal gas law. We capture the effects of variable meniscus curvature and, subsequently, when applicable, movement of triple contact lines. Results are presented for the evolution of the dissolved gas concentration field in the liquid and, when applicable, the time it takes for a meniscus to depin and that for longevity, i.e., the onset of the Cassie to Wenzel state transition. Two configurations are examined; viz., one where an impermeable membrane pressurizes the liquid above the ridges and one where hydrostatic pressure is considered and the top of the liquid is exposed to noncondensible gas.

## 1 Introduction

Superhydrophobic surfaces (SHs) are structured surfaces that can support liquid in the unwetted (Cassie) state. They consist of a hydrophobic substrate on which liquid is supported between protruding structure elements by surface tension and (non-condensible) gas pressure, resulting in both gas and solid contact with the liquid (the vapor pressure of the liquid may also provide support). The distance between the protruding structures can vary but is often of the order of tens to hundreds of microns (see, Refs. [1–3]). A liquid is said to be in the Cassie state when the cavities between the surface elements are unwetted and in the Wenzel state when the liquid fills these cavities. For liquid droplets in the Cassie state, the reduction in total solid–liquid interfacial area results in apparent contact angles approaching 180 deg [4]. In liquid flow in the Cassie state over SHs, this reduction in the solid–liquid interfacial area reduces frictional drag [5]; however, our focus in this paper is on a stationary liquid. Engineering and biological applications of SHs with stationary liquids include antifouling and respiration of underwater insects [6,7], respectively.

For a liquid column that is supported in the Cassie state, the balance of forces along with liquid–vapor interfaces (menisci), assuming they are stationary or move very slowly, is given by the Young–Laplace equation. However, applied pressure and gas diffusion can cause a liquid column that is initially in the Cassie state to transition to the wetted (Wenzel) one [3,4,8–10]. The commonly studied ridge-type SHs with closed cavities are our focus here. Upon submersion in a column of water, the gas in the cavity(s) of SHs with closed cavities is rapidly compressed, in what may be assumed an isothermal or adiabatic process, to a pressure we denote by $pg,c(0+)$ [11,12]. The corresponding local increase in the dissolved gas concentration on the liquid side of the menisci in accordance with Henry's Law causes gas molecules to diffuse into the liquid. The diffusion of gas out of the cavity results in a decrease of the gas pressure as per the ideal gas law. Meniscus deformation can be separated into two regimes. In the first regime, the meniscus must satisfy a force balance as per the Young–Laplace equation and the contact angle (*θ*) increases. As gas continues to diffuse, the contact angle continues to increase and may reach its maximum value, i.e., the advancing contact angle ($\theta a$) [8,11]. We refer to this time as a critical time $(tcr)$. The second regime of meniscus deformation occurs after critical time. In this regime, when additional gas molecules diffuse into the liquid, the volume of gas in a cavity further decreases, and thus menisci depin and slide down the cavity walls [8]. Finally, the menisci may reach the bottom surface (Wenzel state). The duration of the Cassie to Wenzel state transition is known as longevity and is denoted $tf$ [10,12]^{2}.

The upward component of the surface tension force during the metastable Cassie state depends on the shape of the cavity sidewalls [13]. For example, for ridge-type structures that have straight side walls, it remains fixed during movement of the triple contact line. Conversely, for other sidewall profiles, such as cones or reentrant structures, it is dependent on the location of the triple contact line [13]. In some experimental studies, the onset of the Wenzel state has been observed before the menisci fronts are theoretically expected to contact the bottom surface, possibly as a result of imperfections in the cavity walls or the presence of another favorable minimum energy state [8,12,14].

The mechanisms for Cassie to Wenzel state transition can be observed and validated experimentally using, for example, optical [10,14] and acoustic techniques [15]. Lv et al. [8] reported the menisci shape during a pressure- and diffusion-induced transitions on micropore type structures using a confocal microscope. They measured the contact angle and displacement of a triple contact line as a function of time. Their results demonstrate the two regimes of meniscus deformation, initially showing an increase in contact angle followed by sliding into micropores upon achievement of the advancing contact angle. Similarly, Xu et al. [3] used a high-resolution camera to capture the gas diffusion-induced transition on a microscale trench sample. Their method is discussed in detail later. Additionally, several of the experimental studies showed that increasing the degree of gas saturation in the liquid can slow the diffusion process and thus enhance the longevity of SHs [10,12].

Early work on modeling gas diffusion-induced Cassie to Wenzel state transition on ridge-type structures was done by Emami et al. [11]. In their model, the diffusion of dissolved gas into the liquid is proportional to the product of gas pressure in the cavity and an experimentally determined average mass transfer coefficient for air. This model has been used by others to correlate experimental data on gas diffusion-induced Cassie to Wenzel state transition by prescribing an appropriate value for the aforementioned mass transfer coefficient [3,10]. In contrast, for liquid over cylindrical microcavities, Søgaard et al. [14] numerically solved a one-dimensional transient gas diffusion model and successfully correlated their experimental data without the need for experimentally-determined fitting parameters. In their model, the radius of curvature of the meniscus was held constant at the advancing contact angle, ignoring meniscus deformation before critical time. The Young–Laplace equation was used to determine pressure within the cavity and this pressure was coupled to a one-dimensional diffusion equation via Henry's constant.

Previously, we developed a semi-analytical, one-dimensional solution for the evolving dissolved gas concentration field in the liquid, the critical time, and longevity on ridge-type SHs [16]. We captured the effect of a time-dependent boundary condition on the menisci by applying Duhamel's theorem. Meniscus curvature was accounted for in the Young–Laplace equation and ideal gas law. However, in solving the transient diffusion equation, we assumed that the solid fraction of the ridges, defined as the meniscus length divided by ridge pitch, was vanishingly small and the meniscus was flat. In this numerical study, we relax the foregoing assumptions and thus refine our predictions for times required for menisci to depin and transition to the Wenzel state.

Xu et al. [3] experimentally capture the effects we model in this paper. They submerged a single superhydrophobic cavity whose length was much larger than its width in water and imaged the meniscus location over time. Initially, the meniscus was flat but immediately upon submersion the contact angle started to increase as gas diffused into the liquid. In the middle of the long cavity, far from the end walls, the deforming meniscus was unaffected by edge effects and the deformation appeared two-dimensional when viewed from the side. Upon reaching the advancing contact angle, the meniscus depinned from the cavity edges and proceeded to travel into the cavity with a constant radius of curvature. A complimentary model was developed using an empirical parameter to capture meniscus curvature effects. We build a more accurate model than the one implemented in Xu et al. [3] by fully resolving the two-dimensional gas distribution in the liquid, explicitly capturing meniscus curvature and transient, two-dimensional diffusion. Physically our model can be viewed as a two-dimensional slice of a periodic superhydrophobic surface in which the slice is taken in the middle of a long channel such as the one seen in the experiments by Xu et al. [3]. However, we note that Xu et al. [3] focus on a single channel and our model is the first to our knowledge that resolves a fully periodic system of grooves and thus there is no direct experimental data available against which to compare our model. In contrast to our schematic configuration, Xu et al. [3] submerge their single grooves in a cylindrical container. Too, periodic systems of grooves are more common than single grooves in practice due to larger surface coverage.

## 2 Mathematical Model

*Ω*, between the top of a ridge and the (assumed flat) top surface of the liquid,

*s*

_{0}is the ridge height (or cavity depth), 2

*a*is the width of the cavity and 2

*d*is the ridge period such that half of the ridge thickness is (d–a). We additionally define the solid fraction of the ridge as $\varphi =(d\u2212a)/d$. Assuming the pressure on the liquid side of the meniscus is equal to or larger than that on the gas side of it, a stress balance along the meniscus is given by

*g*is the acceleration due to gravity, $ym(x,t)$ is the

*y*-coordinate of the meniscus,

*σ*is the liquid–gas surface tension and

*R*is the radius of curvature of the meniscus. $ym(x,t)$ is a moving boundary and its location must be tracked during our solution. We assume that $h(t)\u226b|ym(x,t)|\u2192h(t)\u2212ym(x,t)\u2248h(t)$; therefore, changes in hydrostatic pressure along the meniscus are negligible and

*R*is independent of

*x*, i.e., the meniscus is a circular arc. Consequently, Eq. (1) simplifies to

^{3}. Also, it follows from the geometry of the problem that

*s*(

*t*) is the distance (a positive quantity) from the top of a ridge to the triple contact line and the ‘+’ and ‘–’ signs preceding the bracketed term correspond to menisci protruding upward and downward, respectively. Henceforth, we restrict our attention to horizontal and downward protruding menisci. As per Fig. 1, the contact angle, i.e., the (positive) angle between the vertical and the tangent to the meniscus at the triple contact line, is between $\pi /2$ and the advancing contact angle. Before critical time, i.e, when $s(t)=0$, the meniscus is pinned and

*R*decreases with time. After the critical time it assumes a constant value corresponding to the advancing contact angle. The geometric relationship between

*R*and

*θ*is

*x*=

*0. The volume of gas per unit depth in one half of a cavity, $Vg,c\u2032(t)$, is given by $Vg,c\u2032(t=0)=a\u2009s0$ plus the integral of $ym(x,t)$ on $0\u2264x\u2264a$. It follows that*

where $mg,c\u2032$ is the residual mass of gas in the cavity per unit depth, $Mg$ is the molecular weight of the gas, $Ru$ is the universal gas constant and *T* is the absolute temperature of the gas.

*t*=

*0, gas pressure in the cavity, $pg,c(0)$, is equal to $p\u221e$ and thus, the meniscus is flat. Figure 1(b) shows a schematic of case study $II$, where we consider the effects of hydrostatic pressure and gas diffusion at the top of the liquid domain. Because of hydrostatic pressure, the initial pressure in the cavity will be higher than that of the atmosphere and the meniscus has an initial radius of curvature of $R(0+)$ [11]. To solve for the initial volume, we assume that the gas pressure in the cavity is initially (i.e., at*

*t*=

*0) $p\u221e$. However, we assume that the gas is instantly compressed to a pressure of $pg,c(0+)$ and a volume per unit depth of $Vg,c\u2032$ due to hydrostatic pressure. Assuming isothermal compression, $pg,c(0+)$ is given by*

The meniscus has been assumed flat $p\u221e$ at the beginning of the isothermal compression process. These are reasonable assumptions when pouring water onto a uniformly-covered textured surface as air will be trapped in the trench at atmospheric pressure and initially the meniscus supports a negligible pressure difference. We note that the meniscus immediately depins if the hydrostatic pressure reduces volume in the cavity below the volume needed to reach the advancing contact angle, i.e., $s(0+)\u22600$ and $tcr=0$, if $Vg,c\u2032(0+)<Vg,c\u2032(\theta =\theta a)$.

*x*=

*0 and*

*x*=

*d*as per Eqs. (9) and (10). The cavity is sealed on all sides; therefore, the gas leaves it only through diffusion across the meniscus. Dissolved gas concentration on the liquid side of the meniscus is determined from Henry's Law as per Eq. (11), where

*H*is Henry's constant. Case study I and case study II only differ in the boundary condition at the top of the domain and share the following boundary conditions.

*s*is the coordinate along with the meniscus. Upon substituting Eq. (19) into Eq. (18) and transforming the surface integral along $ds$ in Eq. (18) to one along $dx$ and integrating over time we obtain

A dimensionless analysis is forgone due to the appearance of numerous dimensionless numbers that appear when scaling the various constraints in this problem.

## 3 Numerical Method

### 3.1 Introduction.

We apply a spectral (Chebyshev collocation) method to numerically solve the transient, two-dimensional diffusion equation, Eq. (8), for the partial density of the dissolved gas in the liquid. This preserves discretization accuracy near curved boundaries. However, once transformed, the diffusion equation and boundary conditions are more complex.

### 3.2 Domain Decomposition and Mapping.

As per Fig. 2, we solve for the dissolved gas concentration field in domain, Ω by decomposing it into an arbitrary number of subdomains, Ω* _{q}*, for $q=1,2,\u2026,Q,Q+1,\u2026,2Q$ for $t\u2264tcr$ and $q=1,2,\u2026,Q,Q+1,\u2026,2Q,2Q+1$ for $t>tcr$ as shown in Fig. 2. The columns of subdomains to the left and right of the triple contact point at $x=a,y=0$ for $t\u2264tcr$ and $x=a,y=\u2212s(t)$ for $t>tcr$ are separated by a (vertical) line emanating from it. There are

*Q*rows and 2 columns of subdomains up until $tcr$, when subdomain $2Q+1$ is added to the bottom of the left column. Straight boundaries of subdomains parallel to the

*x*axis correspond to fixed values of

*y*denoted by

*y*, where $q=0,1,2,\u2026,Q\u22121$, except at the top of the domain, where the value of

_{q}*y*decreases with time and corresponds to $h(t)$. Care must be taken in choosing the height of the topmost subdomains to ensure the moving boundary there can be captured in that subdomain for the entirety of the solution. Heights of the subdomains are chosen so that those closest to the meniscus have aspect ratios near unity. Nodes are aligned on subdomain interfaces. Before critical time, the left, bottom subdomain, $\Omega 1$, is bounded by a deforming boundary, $y=ym(x,t)$, and a fixed one,

*y*=

*y*

_{1}. At critical time, the meniscus depins and the triple contact line moves down the impermeable ridge; therefore, we separate $\Omega 1$ into two subdomains as per Fig. 2(b). Then, subdomain $\Omega 1$ becomes rectangular and bounded from below by

*y*=

*0 and subdomain $\Omega 2Q+1$ contains the rest of the liquid previously in subdomain $\Omega 1$. In general, the subdomains fall into 7 categories based on their applicable boundary conditions and geometry and this is discussed in detail in Appendix A.*

*ξ*and

_{q}*η*are given by the maps

_{q}where $xq,\u2212,\u2009xq,+,\u2009yq,\u2212(x,t)$, and $yq,+(t)$ are the positions of the left, right, lower and upper boundaries, respectively, on subdomain *q* [19]. Here $yq,\u2212$ only depends upon *x* and *t* for subdomain 1 when $t\u2264tcr$ and subdomain $2Q+1$ when $t>tcr$ and $yq,+$ only depends on *t* for subdomains *Q* and 2*Q*.

*x*and

*y*in terms of (ordinary and partial) derivatives of the maps and those of $\rho (\xi ,\eta ,t)$ are

*ξ*,

_{xx}*ξ*,

_{y}*ξ*and

_{yy}*η*equal 0 in all subdomains, the second-order partial derivatives are

_{yy}*η*is only nonzero on subdomain 1 before critical time and subdomain $2Q+1$ after it. Furthermore,

_{xx}*y*of physical subdomains 1,

*Q*and 2

*Q*for $t\u2264tcr$ and

*Q*, 2

*Q*, and $2Q+1$ for $t>tcr$ when

*ξ*and

*η*are held constant [20,21]. It follows that the transient diffusion equation, Eq. (8), in the (Cartesian) $\xi \eta $ plane is

*q*to the dependent variable

*ρ*to emphasize our splitting of our domain into subdomains. Here, $l\u2113(Lq)(\xi q)$ and $ln(Nq)(\eta q)$ are the $\u2113$-th and

*n*-th terms of

*L*-th and

_{q}*N*-th order Lagrange-basis polynomials in

_{q}*ξ*and

_{q}*η*, respectively, as per

_{q}*I*and

_{y}*I*are the $(N+1)\xd7(N+1)$ and $(L+1)\xd7(L+1)$ identity matrices, respectively,

_{x}*D*is an $(L+1)\xd7(L+1)$ Chebyshev differentiation matrix as per

_{x}*D*is the analogous $(N+1)\xd7(N+1)$ Chebyshev differentiation matrix.

_{y}^{4}The symbol ⊗ denotes the tensor product, which produces matrices of the requisite $(L+1)(N+1)\xd7(L+1)(N+1)$ dimensions. The second-order partial derivatives are

### 3.3 Time Discretization.

*k*-th time-step, we march forward with a constant time-step size,

*δt*. For constant

*δt*, $(\rho t)\xi ,\eta $ in Eq. (29) is given by

*k*), and $(k\u22121)$ denote the number of time steps elapsed since

*t*=

*0. The same formula is used to evaluate $dR/dt$ in subdomains 1 and 2*

*Q*+

*1 and $dh/dt$ in subdomains*

*Q*and 2

*Q*. The analogous two-point formula is used for the first time-step. As the contact angle approaches the advancing contact angle before critical time, we reduce the time-step significantly for better precision. This reduces the three-point formula to the form,

where $\delta t\u2032$ is the new time-step size and $dR/dt$ and $dh/dt$ are computed in a similar manner. We note that Eq. (46) simplifies to Eq. (45) for $\delta t\u2032=\delta t$. The first-time-step after critical time was treated with a two-point time discretization before reverting to a 3-point formula.

### 3.4 Transformed Diffusion Equation.

*x*'s and

*y*'s to

*ξ*'s and

*η*'s, respectively, in Eqs. (36), (37), (43), and (44), provides $\rho \xi $, $\rho \eta ,\u2009\rho \xi \xi $, and $\rho \eta \eta $, respectively. It then follows from Eqs. (24)–(27) that

*η*and

*ξ*are $(L+1)(N+1)\xd7(L+1)(N+1)$ diagonal matrices. Expressing the diffusion equation in the form

for the case of three-point time discretization. Here, *I _{q}* is an $(L+1)(N+1)\xd7(L+1)(N+1)$ identity matrix and

*y*has the same dimensions. We note that, in certain subdomains, $\eta x$,

_{t}*η*, $\eta y$, and $yt$ change with time; therefore, at each time-step, we update them using the subdomain-level equations as described in Sec. 3.6.

_{xx}### 3.5 Assembly of Linear System of Equations.

*-Ω*

_{i}*interface. (If Ω*

_{j}*and Ω*

_{i}*do not share an interface, then $A~i,j=0$.) The resulting linear system of equations is*

_{j}In discussing this system of equations we refer to the matrix operator as $Aglobal$, the full solution vector as $\rho global$, and the forcing vector as $Fglobal$. Examples of how boundary nodes are accommodated are as follows. First, to prescribe an (external) homogeneous Neumann condition in *x* on the $ith$ (noncorner) Chebyshev node in the column vector $\rho q$, we set the $ith$ row of $Aq$ equal to the right-side of Eq. (47) and the $ith$ row of $Fq$ equal to 0. Secondly, for the node shared by subdomains 1 and 2 along the vertical line of symmetry on the left-side of the domain (see Fig. 2), we modify rows in $A1,\u2009A2,\u2009F1$, and $F2$ and utilize the relevant rows in $A\u03031,2$ and $A\u03032,1$ to impose continuity of temperature and diffusive flux in the *y* direction. (The choice of the interfacial boundary condition, i.e., flux or concentration matching, assigned to a given subdomain is arbitrary.) We do not impose the homogeneous Neumann condition in *x* at this node. Finally, for a corner node shared by 4 subdomains, we match 3 concentrations and one diffusive flux. Additional information on the subject can be found in the texts by Canuto et al. [22] and Trefethen [24].

### 3.6 Sumdomain Types.

The subdomains fall into 7 groups as follows.

The left, bottom subdomain when $0<t\u2264tcr$ in case study I and $0+<t\u2264tcr$ in case study II bounded by $0<x<a$ and $ym<y<y1$, i.e., $\Omega 1$ as per Fig. 2(a). After critical time, it subdivides into $\Omega 1$ of group (B) and $\Omega 2Q+1$ of group (G) as per Fig. 2(b).

Middle subdomains on $0<x<a$, i.e., from $\Omega 2$ to $\Omega Q\u22121$ and, additionally, $\Omega 1$ for $t>tcr$.

The top subdomain on $0<x<a$, i.e., Ω

._{Q}The bottom subdomain on

*a*<*x*<*d*, i.e., $\Omega Q+1$.Middle subdomains on

*a*<*x*<*d*, from $\Omega Q+2$ to $\Omega 2Q\u22121$.The top subdomain on

*a*<*x*<*d*, i.e., $\Omega 2Q$.The bottom subdomain on $0<x<a$ after critical time, $\Omega 2Q+1$. Before critical time, the portion of the solution vector corresponding to $\Omega 2Q+1$ is set to zero because it's only defined after the triple contact line depins.

In order to formulate the necessary system of equations, we describe boundary conditions and geometric maps, $\xi (x)$ and $\eta (x,y,t)$, for each subdomain group. These are discussed in detail in Appendix A. Additionally, as part of the numerical method, we performed singularity removals on the singular points in the problem. The triple contact point is singular before and after critical time and the point at which the meniscus intersects the cavity wall is singular. The singularity removals were done in order to ensure correct capture of the singularity behavior, as well as increase the speed and accuracy of the solution by removing the need for large mesh densities in the vicinity of the singular points. The singularity removals are discussed in detail in Appendix B.

### 3.7 Solution

#### 3.7.1 Solution Up to Critical Time.

*f*

_{1}is the first term in Eq. (21) (evaluated at $s(t)=s0$) and

*f*

_{2}is the second one. An inner loop, denoted by

*i*is used to update the maps $\eta x,\u2009\eta x$, and $yt$. The first iteration of the inner loop is the solution of the linear system of equations where we assume no motion at the meniscus since the last time-step and is denoted by $\rho 1(tk,1)$. Subsequently, the solution $\rho 1(tk,i=1)$ is then substituted into Eq. (54) to update the pressure of gas in the cavity, as denoted by $pg,c(tk=1,i=2)$. We compute the integral in

*f*

_{2}in Eq. (21) from $tk\u22121$ to

*t*using the Clenshaw–Curtis quadrature method as it preserves accuracy of the Chebyshev spectral method [22]. With this improved estimate of $pg,c(tk,i+1)$, we update the maps and the corresponding terms $\eta x,\u2009\eta y$, and $yt$ in Eq. (52) for use in the next iteration. The boundary condition on the meniscus is updated to $pg,c(tk,i=2)H$. In the next iteration,

_{k}*i*=

*2, $pg,c(tk,i=2)$ is the initial condition to compute a more accurate dissolved gas concentration field in the liquid, $\rho q(tk,i=3)$. Finally, we iterate until a convergence criteria is met*

where *ϵ* is user-defined.

The algorithm is:

Determine the $A\u0303q$ and $F\u0303q$ matrices using solution of the previous iteration as an initial condition.

Numerically solve for $\rho num(tk,i)$ and $\rho sing$ as per Eq. (B9).

Compute $pg,c(tk,i+1)$ as per Eqs. (54).

Update the geometric parameters, $\eta x,\u2009\eta y$, and $yt$. Compute time derivatives using a 3-point formula of the form of Eq. (45).

If the convergence criterion given by Eq. (55) is met, then repeat steps 1 through 4 for the next time step $(k+1)$. If the convergence criterion is not met, proceed to the next iteration (

*i*+ 1) for the current time step (*k*).

*i*=

*1, we approximate the shape of the liquid domain and the boundary condition at the meniscus to be the same as those at $t=0+$*

When critical time has been reached, we fix the gas pressure and allow the meniscus to depin.

#### 3.7.2 After Critical Time.

*t*$(for\u2009\u2009tcr<tk\u2264tcr+\delta t)$, we use an initial guess for the finite displacement of the triple contact line $s(tk,i=1)$. We perform a mass balance on the gas cavity in a similar manner to Eq. (54) to obtain

_{k}For subsequent time steps, we employ the same approach as described in the previous section, except we compute $s(tk,i+1)$ instead of $pg,c(tk,i+1)$.

### 3.8 Longevity.

## 4 Results and Discussion

We implemented the model in MATLAB^{®}. Assuming that the liquid and gas are pure water and nitrogen, respectively, and that liquid rests on a silicon substrate coated with a flouropolymer coating, we set the thermophysical properties to $H=1.82\xd710\u22127\u2009kg/(m3Pa)$ [25], $Dgl=2.01\xd710\u22129m2/s$ [26], $\theta A=110\u2009deg$ [17] and $\sigma =0.073\u2009N/m$ [27] at ambient conditions of $T=25\u2009\xb0C$ and $p\u221e=1\u2009atm$. We note in practice that air would likely be the lgas in the cavity; however, since air is primarily nitrogen we assume the gas is pure nitrogen to restrict attention to a single species, simplifying the model. Additionally, $110\u2009deg$ is the advancing contact angle of water on silicon treated with a flouroploymer [17]. Unless noted, we set the geometric parameters to be $s0=10\u2009\mu m$ and $a=5\u2009\mu m$ with the liquid initially degassed. To generate the computational domain for case study $I$, we partitioned the liquid domain into 4 subdomains (before critical time), i.e., *Q *=* *2, where the interfaces *y _{q}* are located on ${0,3a,h}$. This allocated most of the nodes to the region near the triple point where gradients are largest. For case study II, we partitioned the liquid domain into 6 subdomains (before critical time), i.e.,

*Q*=

*3, where the interfaces*

*y*are located on ${0,3a,3a+3(h0\u22123a)/4,h}$. This allocated nodes close to both the triple point and to the top-most meniscus. Before critical time the time stepping was handled as follows. A time-step of $\delta t=10\u22123\u2009s$ was initially used before critical time. However, as the angle of the deforming meniscus neared the advancing contact angle it was changed to $\delta t=10\u22124\u2009s$ to more accurately capture $tcr$. A time-step of $\delta t=0.01\u2009s$ after critical time was needed to accurately capture $tf$. A mesh dependence study was used to calculate the necessary spatial grid size. The total number of nodes was doubled until $tcr$ varied by less than 0.1 percent. This resulted in

_{q}*L*=

*24 and*

*N*=

*20, a total of 1920 nodes, for case study I and*

*L*=

*36 and*

*N*=

*18, a total 3888 nodes, for case study $II$. For both cases $tf$ varied by less than 0.25 percent between the final two iterations. The difference in grid construction between case study I and case study $II$ was in order to better capture the sharp gradients at the top of the domain that occur in case study $II$.*

### 4.1 Case Study I.

To examine case study $I$ we ran a series of computations that spanned the parameter space of solid fraction (i.e., $\varphi =(d\u2212a)/d$ varied between 0 and 1) over a large range of values of the initial height of the liquid column (*h*_{0}). Figure 3 plots critical time versus solid fraction when *h*_{0}, defined as the initial liquid column height, is fixed at 800 *μm* and the half-width of the gas cavity is *a *=* *5 *μ*m. Clearly, increasing solid fraction dramatically increases gas diffusion out of the cavity because of the increase in the width of the domain (2*d*) leading to lateral diffusion. Of note, critical time is almost entirely independent of chamber height, *h*_{0}, for these specific parameters because the time it takes for gas to diffuse through the entire domain is significantly larger than $tcr$. Also plotted is critical time computed from the quasi-one-dimensional model in our earlier publication [16]. This model is exact in the limit as $\varphi \u21920$ and $\theta a\u219290\u2009deg$ and is far easier to solve, but begins to veer from the exact solution as these parameters change in value. As can be seen in Fig. 4 treating the meniscus as flat as opposed to having an advancing contact angle of $110\u2009deg$ leads to an underestimation of the critical time by 6% as $\varphi \u21920$.

*h*

_{0}for solid fractions ranging from 0.1 to 0.9. Again, it is evident that increasing solid fraction speeds up the diffusive process, this time exhibited by the decrease in $tf$. This plot quantifies the increase in $tf$ with decreasing

*h*

_{0}. The (vertical) dashed lines correspond to the limit at which a full Cassie to Wenzel state transition can occur. They represent the minimum liquid column height that has the capacity to hold $m\u20320\u2212m\u2032f$ mass of gas in a saturated solution, where $m\u2032f$ is the mass in the cavity at time $tf$, i.e.

which is a modification to a result found in our earlier publication, i.e., Kadoko et al. [16]. The computed data is again compared to that of Kadoko et al. [16] and again shows underestimation with final time solutions veering from the 2D solution by about 9%.

Figure 5 plots the volume of gas in the cavity versus time for different solid fractions. The inset shows the results for the critical time and slightly thereafter. Since the cavity width is the same, the effect of different solid fractions can easily be explored. As gas diffuses from the cavity into the liquid, the meniscus deforms and eventually depins, reducing the volume of gas in the cavity. This happens significantly faster for domains with the same cavity width, but larger solid fraction because the total volume for diffusion is much larger, approaching that of a half-space for small $\varphi $. The discontinuity visible in the inset highlights the transition between the two types of meniscus movement: (1) meniscus deformation and (2) sliding of the contact line. It exists because after critical time the pressure in the cavity remains constant and *s* changes, whereas before critical time the pressure varies to accommodate the deforming meniscus. $Vmin$ is the volume of gas in the cavity when the lowest part of the meniscus makes contact with the cavity bottom.

Figure 6 shows contour plots of the dissolved gas concentration field in the lower region of the liquid domain, $\u221210\u2009\mu m<y<15\u2009\mu m$, with $\varphi =0.5$. Figure 6(a) shows *ρ* at $t=0.05\u2009s$, soon after the diffusion process has started, such that the concentration gradient is relatively large and the meniscus has only begun to deform. Figure 6(b) depicts *ρ* at critical time, right before the meniscus begins to slide. Figures 6(c) and 6(d) show *ρ* at $t=25\u2009s$ and $t=tf$, respectively. At this point the concentration gradient has decreased significantly and diffusion rates are low. This helps illuminate why $tf$ is much larger than $tc$.

### 4.2 Case Study II.

Figure 7 plots longevity as a function of cavity depth for a characteristic initial liquid height, $h0=3000\u2009\mu m$ and cavity half-width of $a=5\u2009\mu m$. The first thing to note is that a larger initial column height of $h0=3000\u2009\mu m$ in case study $II$ leads to less longevity, $tf$, than in case study $I$ for $h0=1500\u2009\mu m$ and the same cavity depth, $s0=10\u2009\mu m$, as per Fig. 4. This is because in case study $II$ hydrostatic pressure is taken into account. If the effect of hydrostatic pressure is significant the pressure in the cavity will be higher for the same atmospheric pressure, increasing the diffusive driving force. Secondly, Fig. 7 also shows that increasing cavity depth yields larger critical times and longevity because more gas must diffuse out of the cavity to reach the Wenzel state. This is potential valuable insight for design of these surfaces. Plotted alongside the results is the quasi-one-dimensional solution for case study $II$ from Kadoko et al. [16]. In this case, it fulfills the solution requirement for $\varphi =0$, demonstrating much the same shape as our results.

Figure 8 plots longevity for varying initial gas concentration for a characteristic initial liquid height, $h0=3000\u2009\mu m$. Increasing the initial gas concentration reduces the driving force for diffusion, increasing both critical time and longevity. Longevity grows rapidly with increasing initial concentration. This highlights that initial gas concentration, along with cavity depth and total chamber height, can be tuned to achieve the desired performance of a surface. Again to span the whole parameter space, the results from Kadoko et al. [16] for $thetaa=90\u2009deg$ are included.

## 5 Conclusion

We present a two-dimensional numerical analysis of the gas diffusion-induced Cassie to Wenzel state transition in a quiescent liquid resting on parallel-ridge type SHs. This analysis is representative of the meniscus behavior in the middle of a long three-dimensional meniscus such as experimentally studied by Xu et al. [3]. We demonstrate that the algorithm captures the deforming and moving boundaries of the liquid domain seen in experiments and the time-dependent boundary condition on the meniscus in parallel-ridge type SHs. Too, it computes the dissolved gas concentration field in the liquid. The algorithm can be extended to hierarchical reentrant SHPo surfaces where sidewalls of the structures have nonstraight profiles. Future work needed is resolution of the fully three dimensional problem, capturing wall effects near ends of the long channels, and treating the gas in the cavity as air is such that two species equations need to be resolved. Finally, the lack of experimental data to which to compare our results shows a need for experimental work capturing diffusion effects on periodic surfaces.

## Acknowledgment

This material is based upon work supported in part by the National Science Foundation under Grant No. 1402783 and supported in part by a Google Award in Faculty Research to Professor Marc Hodes. We acknowledge useful discussions with Dr. Simon Game and Dr. Georgios Karamanis.

## Funding Data

Google Award in Faculty Research.

National Science Foundation under (Grant No. 1402783; Funder ID: 10.13039/100000001).

## Footnotes

For hierarchical SHs that are composed of micro-structures fabricated on top of a micro-structured substrate, only the nano-structures are wetted when the menisci reach the bottom surface (partial Wenzel state) [9]. A complete Wenzel state is achieved when the nano-structures are wetted.

Although we ignore the effects of changes in meniscus shape and position on hydrostatic pressure, we do account for them in the rest of our formulation as per below as they are essential to, e.g., calculate the volume of the gas phase.

We have dropped the subscript *q* from *L* and *N* as they are the same in all subdomains.

### Appendix A: Subdomain Types

##### A.1 Subdomain Type (A).

*x*and

*y*are

##### A.2 Subdomain Type (B).

*x*=

*0 (for nodes solely in one subdomain) and continuity of dissolved gas concentration and flux of it on the interfaces share by two subdomains. Thus*

*y*are

*Q*– 1 corner points that are shared by four subdomains each, we prescribe three concentration and one flux continuity condition as per

where $q=1,2,\u2026,Q\u22121$.

##### A.3 Subdomain Type (C).

*, is the topmost one on $0<x<a$. At the top side of it we apply impermeability and constant dissolved gas concentration boundary conditions for case studies I and II, respectively. On the bottom and right sides, we apply dissolved gas concentration and flux continuity, respectively. At the left side we apply symmetry conditions. Thus*

_{Q}*y*with respect to

*t*is

##### A.4 Subdomain Type (D).

*x*is

*y*is

##### A.5 Subdomain Type (E).

*x*=

*d*) and continuity on three internal interfaces. Thus

##### A.6 Subdomain Type (F).

*a*<

*x*<

*d*. The boundary condition at the top of $\Omega 2Q$ is the same as that on Ω

*as per Eqs. (A48) and (A49). On the left and bottom sides, we prescribe dissolved gas concentration and flux continuity as per*

_{Q}##### A.7 Subdomain Type (G).

The partial derivative in *y* with respect to *t* is given Eq. (A10), where $dR/dt=0$, because radius of the meniscus is constant after critical time.

### Appendix B: Singularity Removal

##### B.1 Singularity Removal.

A singularity in diffusive flux exists at the triple contact point (*x *=* a*, *y *=* *0) before critical time. After it, as the triple contact point advances down the ridge, one remains at this location and a second one follows the triple contact point. We remove them from our numerical solution (and, subsequently, reincorporate them into our final one) in the manner of Game et al. [28] to reduce the required mesh density as per the following steps. First, local analyses in the vicinity of the singular points are performed to resolve the forms of the singularities up to undetermined constants coefficients, the singularity strengths. Next, a column is added to $A\u0303global$ that contains the contribution(s) to the diffusion equation and (external) boundary conditions from the solution(s) to the local (singular) problem. This requires that the unknown singularity strength(s) be appended to the solution vector. Finally, a row(s) is added to $A\u0303global$ to bound the derivative of the numerical solution at the singularity(s). The resulting matrix system is one column (and one or two rows) larger than prior to singularity removal, a negligible increase in computational load. Below we discuss the specific implementation of each singularity.

###### B.1.1 Before Critical Time.

*B*is the unknown singularity strength, and the eigencondition reads

*α*is a coefficient that depends on the finite difference scheme used. It assumes the values

for the usual three-point formula, the two-point formula and the three-point formula with nonuniform time-step size, respectively. This is necessary, because $\rho sing$ resolves only the Laplacian terms in the diffusion equation, Eq. (52), and the penultimate and ultimate terms on its left side must thus be added.

*x*=

*0, then*

where *i* indicates the *i*-th index and *y*(*i*) is the *y* location of the *i*-th node. Thus, the solution, $\rho global=\rho num+\rho sing$ satisfies the external boundary conditions.

We remove the singularity from the numerical solution, $\rho num$ via the a row vector $vdiff$ in Eq. (B9). It applies the differential operator $\u2202/\u2202x$ at the (singular) triple contact point in subdomain *Q *+* *1. The end result is a new system of $(2Q+1)(L+1)(N+1)+1$ equations. The element 0 in the right-side vector restricts $\u2202\rho num/\u2202x=0$ and therefore removes the singularity from $\rho num$. The full solution is obtained by solving this new system of equations, constructing both the singular solution and the numerical solution from the solution vector (which includes the singularity strength) and adding them together. Additionally, all integrals and derivatives done in the analysis are split into a singular and nonsingular parts and then summed to ensure accuracy.

###### B.1.2 After Critical Time.

The local problems after critical time, denoted by the subscripts 1 and 2 at $x=a,\u2009y=0$ and at the triple contact point, respectively, again cast in polar ($r,\varphi $) coordinates to facilitate their solution, are shown in Fig. 10. Again, locally, we need only resolve Laplace's equation and decompose our solution into numerically resolved and (singular) analytical components.

*x*=

*a*,

*y*=

*0 most singular in the partial derivative with respect to*

*r*is

*C*is the unknown singularity strength. The singular problem at the triple contact point is

*r*is

where *D* is the unknown singularity strength.

The role of the column vectors $v1,diff$ and $v2,diff$ are analogous to that of $vdiff$ in Eq. (B9) before critical time. Moreover, the row vector $v1,diff$ causes the operator $\u2202/\u2202x$ to act on the element of $\rho num$ corresponding to the bottom left corner of subdomain *Q *+* *1 and the row vector $v2,diff$ causes the operator $\u2202/\u2202y$ to act on the element of $\rho num$ corresponding to the triple contact spot in subdomain $2Q+1$. Both these partial derivatives are set equal to 0. The result is a system of $(2Q+1)(L+1)(N+1)+2$ equations. The solution process is analogous to that before critical time.

#### Nomenclature

##### Greek Symbols

*Α*=coefficient for singularity removal

- $\rho $ =
vector of dissolved gas concentration

*δt*=time step

*ϵ*=convergence criteria

*η*=map to vertical of unit square

- Ω =
domain

- $\varphi $ =
solid fraction, $(d\u2212a)/d$

*ρ*=density ($mol/m2$)

*σ*=surface tension along meniscus ($N/m$)

*θ*=contact angle of penetration into cavity, measured from the horizontal

*ξ*=map to horizontal of unit square

##### Roman

- $A$ =
matrix operator

- $A\u0303i,j$ =
matrix for coupling subdomains

*i*and*j**a*=half-width of the cavity (m)

- $B,\u2009C$ =
singularity strengths

*D*=binary molecular diffusion coefficient ($m2/s$) or differentiation matrix

*d*=half of the ridge period (m)

- $F$ =
forcing vector

*g*=gravitational constant ($m2/s$)

*H*=Henry's constant ($mol/m3\xb7Pa$)

*h*=transient height of the liquid domain (m)

*h*_{0}=initial height of the liquid domain (m)

- $m\u2032$ =
mass per unit depth of gas in cavity ($kg/m2$)

- $Mg$ =
molecular weight of gas ($kg/mol$)

*p*=pressure ($Pa$)

*Q*=number of subdomains in vertical direction

*R*=radius of curvature of meniscus (m)

- $Ru$ =
universal gas constant ($J/mol\xb7K$)

*s*=coordinate along arc spanning meniscus (m)

*s*=transient location (in y) of the triple contact point (m)

*s*_{0}=depth of the cavity (m)

*T*=temperature (K)

*t*=time (s)

- v =
contributions from singularity removal to global matrix

- $V\u2032$ =
volume per unit depth of gas in cavity ($m2$)

*y*=_{m}meniscus location (m)

##### Subscripts

- $\u221e$ =
far-field (atmosphere)

- $a$ =
advancing (contact angle)

- $cr$ =
critical

- $diff$ =
differential operator on the triple contact point

- $f$ =
final

- $g,c$ =
gas in the cavity

- $global$ =
referring to global matrix and vectors

- $gl$ =
gas–liquid

- $g$ =
gas

- $l$ =
liquid

- $num$ =
referring to numerical nonsingular solution

- $q$ =
subdomain number

- $sing$ =
referring to quasi-analytical singular solution

## References

*Brunel University Mathematics Technical Papers Collection*.