## 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. [13]). 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,810]. 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 ($θ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

We consider a 2D quiescent liquid column that is initially in the Cassie state above a ridge-type structured surface as per the two configurations shown in Fig. 1. Our focus is restricted to half of a single cell due to periodicity, as discussed in detail later. According to the saturation table for water and the Kelvin Equation, at 25 deg, in the limiting case of a meniscus with a contact angle of 110 deg, the vapor pressure is about 4% of standard atmospheric pressure. 110 deg is the advancing contact angle of water on silicon treated with a flouroploymer such as PTFE, [17], the scenario considered here. Due to the relatively small effect of vapor pressure, we can ignore the vapor component of the mixture which we refer to as the “gas” phase. Additionally, to simplify the model we consider a single gas species in the cavity, in our case nitrogen. In Fig. 1, $h(t)$ is the height of the portion of the liquid domain, Ω, between the top of a ridge and the (assumed flat) top surface of the liquid, s0 is the ridge height (or cavity depth), 2a is the width of the cavity and 2d is the ridge period such that half of the ridge thickness is (d–a). We additionally define the solid fraction of the ridge as $ϕ=(d−a)/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
$pg,c(t)=p∞+ρlg [h(t)−ym(x,t)]−σR(x,t)$
(1)
where $pg,c$ is the gas pressure in the cavity, $p∞$ is atmospheric pressure, $ρl$ is the density of the liquid, 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)≫|ym(x,t)|→h(t)−ym(x,t)≈h(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
$pg,c(t)=p∞+ρlg h0−σR(t)$
(2)
Fig. 1
Fig. 1
Close modal
where $h0=h(t=0)$3. Also, it follows from the geometry of the problem that
$ym(x,t)=±[R(t)2−x2−R(t)2−a2]−s(t)$
(3)
where 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 $π/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
$R(t)={a cos [π−θ(t)] for 0
(4)
As mentioned before, ridge-type SHs are periodic structures so we can focus our solution on a single half-cell of the surface as shown in Fig. 1. The cavities are closed to the environment and thus the only mechanism for gas transport is diffusion into the liquid. Furthermore, we analyze the liquid domain on the interval $0≤x≤d$ because it's symmetric about x =0. The volume of gas per unit depth in one half of a cavity, $Vg,c′(t)$, is given by $Vg,c′(t=0)=a s0$ plus the integral of $ym(x,t)$ on $0≤x≤a$. It follows that
$Vg,c′(t)=a s0−a s(t)−12{R(t)2 arcsin[aR(t)]−a R(t)2−a2}$
(5)
where $R(t)$ can be represented in terms of gas pressure as per Eq. (2). The gas in the cavity obeys the ideal gas law
$pg,c(t)Vg,c′(t)=mg,c′(t)MgRuT$
(6)

where $mg,c′$ 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.

We consider a scenario where a liquid is slowly, as to not allow any impingement of the surface, poured onto an SH, and then immediately gas begins to diffuse into the liquid. Our initial time is the time immediately after the liquid has been poured. In Fig. 1(a), we depict a case study I, where an impermeable membrane separates the liquid from the pressure of the surrounding atmosphere ($p∞$), assumed to be much greater than the hydrostatic pressure of the liquid column, i.e., $p∞≫ρL g h(t)$. In this case, at t =0, gas pressure in the cavity, $pg,c(0)$, is equal to $p∞$ 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∞$. 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′$ due to hydrostatic pressure. Assuming isothermal compression, $pg,c(0+)$ is given by
$p∞a s0=pg,c(0+)Vg,c′(0+)$
(7)

The meniscus has been assumed flat $p∞$ 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+)≠0$ and $tcr=0$, if $Vg,c′(0+).

We neglect the advective mass transport of the dissolved gas species, a valid assumption in analogous mass transfer problems [18]. Thus, the two-dimensional, transient diffusion equation governs the dissolved gas concentration field as per
$∂ρg∂t=Dgl(∂2ρg∂x2+∂2ρg∂y2)$
(8)
where $ρg(x,y,t)$ is the partial density of the dissolved gas in the liquid and $Dgl$ is the binary molecular diffusivity of the dissolved gas in the liquid. We prescribe symmetry boundary conditions at 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.
$∂∂xρg(0,ym(x,t)0)=0$
(9)
$∂∂xρg(d,00)=0$
(10)
$ρg(0
(11)
$∂∂yρg(a0)=0$
(12)
$∂∂xρg(a,−s(t)
(13)
For case study $I$, we assume that the top of the domain is impermeable to gas as per
$∂∂yρg(00)=0, case study I$
(14)
For case study $II$, we assume that the top of the domain is open to the environment such that the dissolved gas concentration on the liquid side of the top of the domain given by Henry's Law as per
$ρg(00)=p∞H, case study II$
(15)
The initial condition is uniformly dissolved gas concentration as per
$ρg(x,y,0)=ρg,0$
(16)
Recall that, for case studies I and II, the initial pressures of gas in the cavity are $p∞$ and $pg,c(0+)$, respectively. In the liquid domain, the total mass of the solvent species is preserved. Moreover, since the partial molar volume of the dissolved gas is much smaller than that of the liquid, the total change in volume of the liquid domain is negligible. It follows from the geometry of the problem that the height of the liquid above the ridge is
$h(t)=h0−ads(t)−12d{R(t)2 arcsin[aR(t)]−a R(t)2−a2}$
(17)
We bound the gas in the cavity by a deforming or moving control volume (of unit depth), $cvg(t)$. The only flux of gas through the corresponding deforming control surface, $csg(t)$, is due to the diffusion of it along with the meniscus. Applying conservation of mass to the gas in the cavity yields
$0=ddt∫cvg(t)ρg,c dA−Dgl∫csg(t)∇ρg·n̂ ds$
(18)
where $ρg,c=mg,c′/Vg,c′$ is the density of gas in the cavity as per Eq. (6), assumed uniform because of the high molecular diffusivity of gases, $dA$ is a differential area in the domain, $n̂$ is the unit vector normal to the meniscus and pointing into the liquid as per
$n̂=−∂ym/∂x i+j1+(∂ym/∂x)2$
(19)
and 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
$0=mg,c′(t)−mg,c′(0)−Dgl∫0t∫0a(∂ρg∂y|y=ym−∂ym∂x∂ρg∂x|y=ym) dx dτ$
(20)
Here, $mg,c′(t)$ and $mg,c′(t=0)$ may be expressed as functions of the pressure and volume of gas in the cavity as per Eq. (6), the ideal gas law. Moreover, the volume of gas in the cavity, $Vg,c′(t)$, is a function of its pressure, from which the radius of curvature of the meniscus follows from Eq. (2), and the displacement of the triple contact line. It follows that Eq. (20) is equivalent to
$0=MgRuT{pg,c(t)Vg′[pg,c(t),s(t)]−p∞s0 a}−Dgl∫0t∫0a[∂ρg∂y−∂ym∂x∂ρg∂x]y=ym dx dτ$
(21)

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,…,Q,Q+1,…,2Q$ for $t≤tcr$ and $q=1,2,…,Q,Q+1,…,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≤tcr$ and $x=a,y=−s(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 yq, where $q=0,1,2,…,Q−1$, except at the top of the domain, where the value of 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, $Ω1$, is bounded by a deforming boundary, $y=ym(x,t)$, and a fixed one, y = y1. At critical time, the meniscus depins and the triple contact line moves down the impermeable ridge; therefore, we separate $Ω1$ into two subdomains as per Fig. 2(b). Then, subdomain $Ω1$ becomes rectangular and bounded from below by y =0 and subdomain $Ω2Q+1$ contains the rest of the liquid previously in subdomain $Ω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.

Fig. 2
Fig. 2
Close modal
We map each subdomain to the canonical $[−1,1]×[−1,1]$ (square) domain suitable for Chebyshev collocation and obtain machinery to effectively accommodate meniscus curvature. The transformed variables ξq and ηq are given by the maps
$ξq(x)=2x−xq,−xq,+−xq,−−1$
(22)
$ηq(x,y,t)=2y−yq,−(x,t)yq,+(t)−yq,−(x,t)−1$
(23)

where $xq,−, xq,+, yq,−(x,t)$, and $yq,+(t)$ are the positions of the left, right, lower and upper boundaries, respectively, on subdomain q [19]. Here $yq,−$ only depends upon x and t for subdomain 1 when $t≤tcr$ and subdomain $2Q+1$ when $t>tcr$ and $yq,+$ only depends on t for subdomains Q and 2Q.

Next, we transform the dependent variable from $ρ(x,y,t)$ to $ρ[ξ(x),η(x,y,t),t]$. The first-order partial derivatives of $ρ(x,y,t)$ with respect to x and y in terms of (ordinary and partial) derivatives of the maps and those of $ρ(ξ,η,t)$ are
$ρx=ξxρξ+ηxρη$
(24)
$ρy=ηyρη$
(25)
where the various derivatives are denoted by subscripts. Noting that ξxx, ξy, ξyy and ηyy equal 0 in all subdomains, the second-order partial derivatives are
$ρxx=ξx2ρξξ+2ξxηxρξη+ηx2ρηη+ηxxρη$
(26)
$ρyy=ηy2ρηη$
(27)
Too, ηxx is only nonzero on subdomain 1 before critical time and subdomain $2Q+1$ after it. Furthermore,
$ρt|ξ,η=ρt|x,y+ρηηyyt$
(28)
where the second term on the right side accounts for the deformation rates in y of physical subdomains 1, Q and 2Q for $t≤tcr$ and Q, 2Q, 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) $ξη$ plane is
$ρt|ξ,η−ρηηyyt=Dgl[ξx2ρξξ+2ξxηxρξη+(ηx2+ηy2)ρηη+ηxxρη]$
(29)
The dissolved gas concentration profile in each subdomain is approximated as a Lagrangian polynomial as
$ρq(ξq,ηq)=∑ℓ=0Lq ∑n=0Nqlℓ(Lq)(ξq) ln(Nq)(ηq)ρq,ℓ,n$
(30)
where we have added the subscript q to the dependent variable ρ to emphasize our splitting of our domain into subdomains. Here, $lℓ(Lq)(ξq)$ and $ln(Nq)(ηq)$ are the $ℓ$-th and n-th terms of Lq-th and Nq-th order Lagrange-basis polynomials in ξq and ηq, respectively, as per
$lℓ(Lq)(ξq)=∏k=0,k≠ℓLqξq−ξq.kξq,ℓ−ξq,k for k=0,1,2,…,Lq$
(31)
$ln(Nq)(η)=∏k=0,k≠nNqηq−ηq,kηq.n−ηq,k for k=0,1,2,…,Nq$
(32)
and $ρq.l,n$ are unknown values of the dissolved gas concentration at predetermined interpolation points, $ξq,ℓ$ and $ηq.n$. These are chosen as the Chebyshev points as per
$ξq,ℓ=−cos(πLℓ), for ℓ=0,1,2,…,Lq$
(33)
$ηq,n=−cos(πNn), for n=0,1,2,…,Nq$
(34)
such that they vary from –1 to 1 and they are collocated with points in the physical domain as per the inverses of the maps given by Eqs. (22) and (23). The uneven spacings of the Chebyshev points reduce Runge's phenomenon at the edges of the computational domain [22]. We note that $lℓ(Lq)(ξq)$ and $ln(Nq)(ηq)$ equal 1 where $ξq=ξq,ℓ$ and $ηq=ηq,n$, respectively, and 0 on all other Chebyshev points. The dissolved gas concentrations at the Chebyshev nodes in each subdomain are expressed in the form of a $(Lq+1)(Nq+1)$ column vector as per
$ρq=[ρ0,0,…,ρLq,0,ρ0,1,…,ρLq,1,…,ρ0,Nq,…,ρLq,Nq]T$
(35)
Given Eq. (30), the column vectors of first-order partial derivatives of the dissolved gas concentration field are [23]
$∂ρq∂x=[Iy⊗Dx]ρq$
(36)
$∂ρq∂y=[Dy⊗Ix]ρq$
(37)
Here, Iy and Ix are the $(N+1)×(N+1)$ and $(L+1)×(L+1)$ identity matrices, respectively, Dx is an $(L+1)×(L+1)$ Chebyshev differentiation matrix as per
$Di,j=c¯ic¯j(−1)i+jxi−xj, i≠j$
(38)
$Di,i=−xi2(1−xi2), 1≤i≤L−1$
(39)
$D0,0=−2L2+16$
(40)
$DL,L=2L2+16$
(41)
where the coefficients $c¯i$ and $c¯j$ are
$c¯ℓ={2 ℓ=0 or ℓ=L1 1≤ℓ≤L−1$
(42)
and Dy is the analogous $(N+1)×(N+1)$ Chebyshev differentiation matrix.4 The symbol ⊗ denotes the tensor product, which produces matrices of the requisite $(L+1)(N+1)×(L+1)(N+1)$ dimensions. The second-order partial derivatives are
$∂2ρq∂x2=[Iy⊗Dx]2ρq$
(43)
$∂2ρq∂y2=[Dy⊗Ix]2ρq$
(44)

### 3.3 Time Discretization.

We discretized the rate of change of dissolved gas concentration at the Chebyshev nodes using a 3-point formula. Then, we compute the dissolved gas concentration using an implicit scheme. On a given k-th time-step, we march forward with a constant time-step size, δt. For constant δt, $(ρt)ξ,η$ in Eq. (29) is given by
$[∂ρ(k+1)∂t]ξ,η=3ρ(k+1)−4ρ(k)+ρ(k−1)2δt+O(δt2)$
(45)
where the superscripts $(k+1)$, (k), and $(k−1)$ denote the number of time steps elapsed since t =0. The same formula is used to evaluate $dR/dt$ in subdomains 1 and 2Q +1 and $dh/dt$ in subdomains Q and 2Q. 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,
$∂ρ(k+1)∂t=[δt(δtδt′+1)]−1{(2+δtδt′)δtδt′ρ(k+1)−(1+δtδt′)2ρ(k)+ρ(k−1)}+O[(δt)2]$
(46)

where $δt′$ 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 $δt′=δ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.

The transformed diffusion equation is expressed in matrix form for each subdomain. Changing the x's and y's to ξ's and η's, respectively, in Eqs. (36), (37), (43), and (44), provides $ρξ$, $ρη, ρξξ$, and $ρηη$, respectively. It then follows from Eqs. (24)(27) that
$∂ρ(k+1)∂x=[ξx(Iη⊗Dξ)+ηx(Dη⊗Iξ)]ρ(k+1)$
(47)
$∂ρ(k+1)∂y=[ηy(Dη⊗Iξ)]ρ(k+1)$
(48)
$∂2ρ(k+1)∂x2={[ξx(In⊗Dξ)+ηx(Dη⊗Iξ)]2+ηxx(Dη⊗Iξ)}ρ(k+1)$
(49)
$∂2ρ(k+1)∂y2=[ηy(Dη⊗Iξ)]2ρ(k+1)$
(50)
where the superscript $(k+1)$ denotes the dissolved gas concentration in the liquid at time $t=(k+1)δt$ and the various partial derivatives of η and ξ are $(L+1)(N+1)×(L+1)(N+1)$ diagonal matrices. Expressing the diffusion equation in the form
$Aqρq(k+1)=Fq$
(51)
where $Aq$ and $Fq$ are an $(L+1)(N+1)×(L+1)(N+1)$ matrix and an $(L+1)(N+1)×1$ vector, respectively, it becomes
$(Dgl{[ξx(Iη⊗Dξ)+ηx(Dη⊗Iξ)]2+ηxx(Dη⊗Iξ)+[ηy(Dη⊗Iξ)]2}+ηyyt(Dη⊗Iξ)−32δtIq)ρ(k+1)=ρ(k−1)−4ρ(k)2δt$
(52)

for the case of three-point time discretization. Here, Iq is an $(L+1)(N+1)×(L+1)(N+1)$ identity matrix and yt has the same dimensions. We note that, in certain subdomains, $ηx$, ηxx, $η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.

### 3.5 Assembly of Linear System of Equations.

Rows in matrices $A1−A2Q+1$ and the corresponding rows in $F1−F2Q+1$ corresponding to internal nodes are assembled according to Eq. (52). We apply the necessary boundary conditions by modifying the $2(Lq+Nq)$ rows in them corresponding to external boundaries or internal ones where subdomains intersect and, denoting the revised matrices by $Ã1−Ã2Q+1$, and by generating the matrices $Ãi,j$ and $Ãj,i$. The latter are sparse matrices necessary to impose the continuity of dissolved gas concentration or flux at the Chebyshev nodes located on an Ωij interface. (If Ωi and Ωj do not share an interface, then $A~i,j=0$.) The resulting linear system of equations is
$[Ã1⋯Ã1,q⋯Ã1,2Q+1⋮⋱⋮Ãq,1ÃqÃq,2Q+1⋮⋱⋮Ã2Q+1,1⋯Ã2Q+1,q⋯Ã2Q+1][ρ1(k+1)⋮ρq(k+1)⋮ρ2Q+1(k+1)]=[F̃1⋮F̃q⋮F̃2Q+1]$
(53)

In discussing this system of equations we refer to the matrix operator as $Aglobal$, the full solution vector as $ρ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 $ρ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, A2, F1$, and $F2$ and utilize the relevant rows in $Ã1,2$ and $Ã2,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.

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

2. Middle subdomains on $0, i.e., from $Ω2$ to $ΩQ−1$ and, additionally, $Ω1$ for $t>tcr$.

3. The top subdomain on $0, i.e., ΩQ.

4. The bottom subdomain on a < x < d, i.e., $ΩQ+1$.

5. Middle subdomains on a < x < d, from $ΩQ+2$ to $Ω2Q−1$.

6. The top subdomain on a < x < d, i.e., $Ω2Q$.

7. The bottom subdomain on $0 after critical time, $Ω2Q+1$. Before critical time, the portion of the solution vector corresponding to $Ω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, $ξ(x)$ and $η(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.

We iteratively solve for the pressure of the gas in the cavity and the dissolved gas concentration in the liquid. At a given $kth$ time-step we numerically solve for an improved estimate of the pressure of gas in the cavity, $pg,c(tk,i+1)$, as per
$0=f1[pg,c(tk,i+1)]+f2[pg,c(tk,i),ρ1(tk,i)]$
(54)
where f1 is the first term in Eq. (21) (evaluated at $s(t)=s0$) and f2 is the second one. An inner loop, denoted by i is used to update the maps $ηx, η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 $ρ1(tk,1)$. Subsequently, the solution $ρ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 f2 in Eq. (21) from $tk−1$ to tk 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 $ηx, η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, i =2, $pg,c(tk,i=2)$ is the initial condition to compute a more accurate dissolved gas concentration field in the liquid, $ρq(tk,i=3)$. Finally, we iterate until a convergence criteria is met
$|pg,c(tk,i+1)−pg,c(tk,i)pg,c(tk,i+1)|≤ϵ$
(55)

where ϵ is user-defined.

The algorithm is:

1. Determine the $Ãq$ and $F̃q$ matrices using solution of the previous iteration as an initial condition.

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

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

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

5. 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).

At $t=0+$, the shape of the domain and dissolved gas concentration in the liquid are given by the initial conditions. At the end of the first time-step and first (inner) iteration, $t=δt$ and 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+$
$pg,c(tk=1,i=1) ≈ {p∞,for case study Ipg,c(0+) ,for case study II$
(56)

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

#### 3.7.2 After Critical Time.

After critical time, the gas pressure and radius of curvature of the meniscus are constant; therefore, $pg,c(t≥tcr)=pcr, θ(t≥tcr)=θA$ and $R(t≥tcr)=Rcr$. However, the meniscus depins such that $s(t)$ is no longer constant and drives the variation of $ym(x,t)$ in the $y2Q+1$ to $η2Q+1$ map. At the first time-step after critical time, tk$(for tcr, 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
$0=f1[s(tk,i+1)]+f2[s(tk,i),ρ1(tk,i)]$
(57)

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.

We assume that the transition to the Wenzel state occurs when the meniscus reaches the bottom of the gas cavity as depicted in Fig. 1, i.e., $ym(0,tf)=−s0$, where $tf$ is the final time or longevity. In a limiting case of $ϕs→0$ and $θA→90 deg$, the solution to $tf$ approaches that of 1D model described in Kadoko et al. [16]. Therefore, in the 1D model, mass conservation of the gas in the cavity yields
$0=p∞s0MgRuT−ρLgH·Dgltf+2hπ2∑n=1∞(1n2{ρLg hH+(p∞H−ρg, 0)[1−cos(πn)]}[exp(−βn2Dgltf)−1])$
(58)

## 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×10−7 kg/(m3Pa)$ [25], $Dgl=2.01×10−9m2/s$ [26], $θA=110 deg$ [17] and $σ=0.073 N/m$ [27] at ambient conditions of $T=25 °C$ and $p∞=1 atm$. 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 deg$ 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 μm$ and $a=5 μ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 yq 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 yq are located on ${0,3a,3a+3(h0−3a)/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 $δt=10−3 s$ was initially used before critical time. However, as the angle of the deforming meniscus neared the advancing contact angle it was changed to $δt=10−4 s$ to more accurately capture $tcr$. A time-step of $δt=0.01 s$ 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 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., $ϕ=(d−a)/d$ varied between 0 and 1) over a large range of values of the initial height of the liquid column (h0). Figure 3 plots critical time versus solid fraction when h0, 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 (2d) leading to lateral diffusion. Of note, critical time is almost entirely independent of chamber height, h0, 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 $ϕ→0$ and $θa→90 deg$ 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 deg$ leads to an underestimation of the critical time by 6% as $ϕ→0$.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Figure 4 plots longevity ($tf$) versus h0 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 h0. 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′0−m′f$ mass of gas in a saturated solution, where $m′f$ is the mass in the cavity at time $tf$, i.e.
$hmin=(1−ϕ)(m′0−m′f)a(p∞H−ρg,0)$
(59)

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 $ϕ$. 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.

Fig. 5
Fig. 5
Close modal

Figure 6 shows contour plots of the dissolved gas concentration field in the lower region of the liquid domain, $−10 μm, with $ϕ=0.5$. Figure 6(a) shows ρ at $t=0.05 s$, 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 s$ 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$.

Fig. 6
Fig. 6
Close modal

### 4.2 Case Study II.

Figure 7 plots longevity as a function of cavity depth for a characteristic initial liquid height, $h0=3000 μm$ and cavity half-width of $a=5 μm$. The first thing to note is that a larger initial column height of $h0=3000 μm$ in case study $II$ leads to less longevity, $tf$, than in case study $I$ for $h0=1500 μm$ and the same cavity depth, $s0=10 μ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 $ϕ=0$, demonstrating much the same shape as our results.

Fig. 7
Fig. 7
Close modal

Figure 8 plots longevity for varying initial gas concentration for a characteristic initial liquid height, $h0=3000 μ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 deg$ are included.

Fig. 8
Fig. 8
Close modal

## 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

2

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.

3

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.

4

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).
The impermeability boundary conditions and Henry's Law, i.e., Eqs. (9) and (11), respectively, apply on exterior faces and on the $Ω1$-$ΩQ+1$ and $Ω1$-$Ω2$ interfaces, we match fluxes. Thus
$∂∂xρ1(0,y,t)=0$
(A1)
$ρ1[x,ym(x,t)]=pg,c(t) H$
(A2)
$∂∂xρ1(a,y,t)=∂∂xρQ+1(a,y,t)$
(A3)
$∂∂yρ1(x,y1,t)=∂∂yρ2(x,y1,t)$
(A4)
We transform these conditions to be in terms of ξ and η as per Eqs. (24) and (25) using the maps
$ξ q(x)=2xa−1, for 0
(A5)
$η1(x,y,t)=2y−ym(x,t)y1−ym(x,t)−1, for t
(A6)
where $ξ q(x)$ is valid for groups (A), (B), (C), and (G) and $η1(x,y,t)$ is only valid for the interval $0 because subdomain $Ω1$ splits into an $a×y1$ rectangle (also denoted by $Ω1$ and belonging to group (B)) and $Ω2Q+1$ after it. The derivative of $ξ q(x)$ is given by
$dξ qdx =2a, for 0
(A7)
The partial derivatives of $η1(x,y,t)$ with respect to x and y are
$∂η1∂x=η−1y1−ym(x,t)xR(t)2−x2$
(A8)
$∂η1∂y=2y1−ym(x,t)$
(A9)
Also, required, for groups (A) and (G), is the partial derivative of y with respect to time, which follows from Eqs. (3) and (23) as
$(∂y(A,G)∂t)ξ,η=1−η2{[R(t)R(t)2−a2−R(t)R(t)2−x2]dRdt−dsdt}, for 0
(A10)
##### A.2 Subdomain Type (B).
The boundary conditions for subdomains in group (B) consist of symmetry on 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
$∂∂xρq(0,y,t)=0$
(A11)
$∂∂xρq(a,y,t)=∂∂xρq+Q(a,y,t)$
(A12)
$ρq(x,yq−1,t)=ρq−1(x,yq−1,t)$
(A13)
$∂∂yρq(x,yq,t)=∂∂yρq+1(x,yq,t)$
(A14)
where
$q={2,3,…,Q−1, for 0
(A15)
The vertical map and its derivative with respect to y are
$η q(y)=2y−yq−1yq−yq−1−1$
(A16)
$∂η q∂y=2yq−yq−1$
(A17)
As for the Q – 1 corner points that are shared by four subdomains each, we prescribe three concentration and one flux continuity condition as per
$ρq(a,yq,t)=ρq+1(a,yq,t)$
(A18)
$ρq+1(a,yq,t)=ρq+Q+1(a,yq,t)$
(A19)
$ρq+Q+1(a,yq,t)=ρq+Q(a,yq,t)$
(A20)
$∂∂xρq+Q(a,yq,t)=∂∂xρq(a,yq,t)$
(A21)

where $q=1,2,…,Q−1$.

##### A.3 Subdomain Type (C).
Subdomain type (C), i.e., ΩQ, is the topmost one on $0. 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
$∂∂xρQ(0,y,t)=0$
(A22)
$∂∂xρQ(a,y,t)=∂∂xρq+Q(a,y,t)$
(A23)
$ρQ(x,yQ−1,t)=ρQ−1(x,yQ−1,t)$
(A24)
$∂∂yρQ(x,h(t),t)=0, for case study I$
(A25)
$ρQ(x,h(t),t)=p∞H, for case study II$
(A26)
The vertical map and its derivative are
$ηQ(y,t)=2y−yQ−1h(t)−yQ−1−1$
(A27)
$∂ηQ∂y=2h(t)−yQ−1$
(A28)
For groups (C) and (F), the rate of change of y with respect to t is
$(∂y(C,F)∂t)ξ,η=12(1+η)dhdt$
(A29)
##### A.4 Subdomain Type (D).
Subdomain type (D), i.e., $ΩQ+1$, is directly above the ridge. The symmetry and impermeability boundary conditions apply to the exterior. We apply dissolved gas concentration and flux continuity on the left and top interfaces. Thus
$ρQ+1(a,y,t)=ρ1(a,y,t)$
(A30)
$∂∂xρQ+1(d,y,t)=0$
(A31)
$∂∂yρQ+1(x,0,t)=0$
(A32)
$∂∂yρQ+1(x,y1,t)=∂∂yρQ+2(x,y1,t)$
(A33)
The horizontal and vertical maps are
$ξ q(x)=2x−ad−a−1, for a
(A34)
$ηQ+1(y,t)=2yy1−1$
(A35)
where Eq. (A34) is valid for subdomain groups (D), (E) and (F). The derivative of $ξ q$ with respect to x is
$∂ξq∂x=2d−a$
(A36)
and the derivative of the vertical map with respect to y is
$∂ηQ+1∂y=2y1$
(A37)
##### A.5 Subdomain Type (E).
Similar to subdomain group (B), boundary conditions for group (E) consist of a symmetry condition on one exterior side (x = d) and continuity on three internal interfaces. Thus
$ρq(a,y,t)=ρq−Q(a,y,t)$
(A38)
$∂∂xρq(d,y,t)=0$
(A39)
$ρq(x,yq−Q−1,t)=ρq−1(x,yq−Q−1,t)$
(A40)
$∂∂yρq(x,yq−Q,t)=∂∂yρq+1(x,yq−Q,t)$
(A41)
where
$q=Q+2,Q+3,…,2Q−1$
(A42)
The map used for subdomains of type (E) is
$η q(y)=2y−yq−Q−1yq−Q−yq−Q−1−1$
(A43)
where
$∂η q(y)∂y=2yq−Q−yq−Q−1$
(A44)
##### A.6 Subdomain Type (F).
Subdomain group (F) is the topmost portion of the liquid domain on a < x < d. The boundary condition at the top of $Ω2Q$ is the same as that on ΩQ as per Eqs. (A48) and (A49). On the left and bottom sides, we prescribe dissolved gas concentration and flux continuity as per
$ρ2Q(a,y,t)=ρQ(a,y,t)$
(A45)
$∂∂xρ2Q(d,y,t)=0$
(A46)
$∂∂yρ2Q(x,yQ−1,t)=∂∂yρ2Q−1(x,yQ−1,t)$
(A47)
$∂∂yρ2Q(x,h(t),t)=0, for case study I$
(A48)
$ρ2Q(x,h(t),t)=p∞H, for case study II$
(A49)

The map in the vertical direction and its derivatives are the same as those given in Eqs. (A27)(A28).

##### A.7 Subdomain Type (G).
Subdomain type (G), i.e., $Ω2Q+1$ exists after critical time as shown in Fig. 2(b). The boundary conditions are
$∂∂xρ2Q+1(0,y,t)=0$
(A50)
$∂∂xρ2Q+1(a,y,t)=0$
(A51)
$ρ2Q+1(x,ym(x,t),t)=pcr H$
(A52)
$∂∂yρ2Q+1(x,0,t)=∂∂yρ1(x,0,t)$
(A53)
The vertical mapping function is
$η2Q+1(x,y,t)=1−2yym(x,t)$
(A54)
The first-order partial derivatives of $η2Q+1(x,y,t)$ are
$∂η2Q+1∂x=1−ηym(x,t)xR(t)2−x2$
(A55)
$∂η2Q+1∂y=−2ym(x,t)$
(A56)

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̃global$ 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̃global$ 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.
The local problem in the vicinity of the triple contact point, cast in polar ($r,ϕ$) coordinates to facilitate its solution, is shown in Fig. 9. Only the Laplacian terms in the diffusion equation are singular; therefore, we need only resolve Laplace's equation in our local analysis. We decompose our solution into a numerically-resolved component, $ρnum$, and an analytically-resolved singular one as per
$ρglobal=ρnum+ρsing$
(B1)
Fig. 9
Fig. 9
Close modal
We impose the in homogemenous boundary condition along the meniscus as per $ρnum=Hpg,c$. Then, both boundary conditions on the local problem are homogeneous and it reads
$∇2ρsing=0$
(B2)
$∂ρsing∂ϕ|ϕ=0=0$
(B3)
$ρsing|ϕ=θ+π/2=0.$
(B4)
The solution is
$Brλ cos(λϕ)$
(B5)
where B is the unknown singularity strength, and the eigencondition reads
$cos [λ(θ+π/2)]=0.$
(B6)
We need only consider the strongest singularity such that
$ρsing=Brλ0 cos(λ0ϕ)$
(B7)
where $λ0=π/(2θ+π)$. In Cartesian coordinates this result is
$ρsing=B(x2+y2)1/(2λ0) cos [λ0 arctan(y/x)]$
(B8)
We subtract $ρsing$ from the numerical problem by recasting the overall problem as
$[Aglobalvsingvdiff0][ρnumB]=[Fglobal0]$
(B9)
Here, on an interior point, $vsing,i$ takes the form
$vsing,i=yt×1B∂ρsing∂y|x=x(i),y=y(i)−α1Bρsing|x=x(i),y=y(i)$
(B10)
where α is a coefficient that depends on the finite difference scheme used. It assumes the values
$α={32δt, 1δt, [δt(δtδt′+1)](2+δtδt′)δtδt′}$
(B11)

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 $ρ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.

When matching flux or concentration on adjacent subdomains singularity subtraction is not necessary and thus the appropriate element of $vsing$ is set to 0. Finally, for points on external boundaries, $vsing$ is populated with the appropriate contribution of $ρsing$. For example, when the $ith$ row of $Aglobal$ corresponds to a node operated by the symmetry boundary condition at x =0, then
$vsing,i=1B∂ρsing∂x|x=0,y=y(i)$

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

We remove the singularity from the numerical solution, $ρnum$ via the a row vector $vdiff$ in Eq. (B9). It applies the differential operator $∂/∂x$ 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 $∂ρnum/∂x=0$ and therefore removes the singularity from $ρ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, y=0$ and at the triple contact point, respectively, again cast in polar ($r,ϕ$) 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.

Fig. 10
Fig. 10
Close modal
A separation of variables solution shows that the form of the solution at x = a, y =0 most singular in the partial derivative with respect to r is
$ρsing,1=Cr2/3 cos(ν1ϕ)$
(B12)
where C is the unknown singularity strength. The singular problem at the triple contact point is
$∇2ρsing,2=0$
(B13)
Subject to
$∂ρsing,2∂ϕ|ϕ2=π/2=0$
(B14)
$ρsing,2|ϕ2=θa+π/2=0.$
(B15)
It follows from the separation of variables and trigonometric identities that the form of the solution at the triple contact line that is most singular in the partial derivative with respect to r is
$ρsing,2=r2π/(2θa)[D sin(πϕ22θa)+Ecot(π2θa)cos(πϕ22θa)]$
(B16)

where D is the unknown singularity strength.

Subtracting $ρsing,1$ and $ρsing,2$ from the numerical problem changes the overall one to
$[Aglobalv1,singv2,singv1,diff00v2,diff00][ρnumCD]=[Fglobal00.]$
(B17)

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 $∂/∂x$ to act on the element of $ρnum$ corresponding to the bottom left corner of subdomain Q +1 and the row vector $v2,diff$ causes the operator $∂/∂y$ to act on the element of $ρ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
Greek Symbols
Α =

coefficient for singularity removal

$ρ$ =

vector of dissolved gas concentration

δt =

time step

ϵ =

convergence criteria

η =

map to vertical of unit square

Ω =

domain

$ϕ$ =

solid fraction, $(d−a)/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
Roman
$A$ =

matrix operator

$Ãi,j$ =

matrix for coupling subdomains i and j

a =

half-width of the cavity (m)

$B, C$ =

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·Pa$)

h =

transient height of the liquid domain (m)

h0 =

initial height of the liquid domain (m)

$m′$ =

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·K$)

s =

coordinate along arc spanning meniscus (m)

s =

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

s0 =

depth of the cavity (m)

T =

temperature (K)

t =

time (s)

v =

contributions from singularity removal to global matrix

$V′$ =

volume per unit depth of gas in cavity ($m2$)

ym =

meniscus location (m)

##### Subscripts
Subscripts
$∞$ =

far-field (atmosphere)

$a$ =

$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

1.
Lee
,
C.
, and
Kim
,
C. J.
,
2009
, “
Maximizing the Giant Liquid Slip on Superhydrophobic Microstructures by Nanostructuring Their Sidewalls
,”
Langmuir
,
25
(
21
), pp.
12812
12818
.10.1021/la901824d
2.
Bolognesi
,
G.
,
Cottin-Bizonne
,
C.
, and
Pirat
,
C.
,
2014
, “
Evidence of Slippage Breakdown for a Superhydrophobic Microchannel
,”
Phys. Fluids
,
26
(
8
), p.
082004
.10.1063/1.4892082
3.
Xu
,
M. C.
,
Sun
,
G. U.
, and
Kim
,
C. J.
,
2014
, “
Infinite Lifetime of Underwater Superhydrophobic States
,”
Phys. Rev. Lett.
,
113
(
13
), p.
136103
.10.1103/PhysRevLett.113.136103
4.
Lafuma
,
A.
, and
Quere
,
D.
,
2003
, “
Superhydrophobic States
,”
Nat. Mater.
,
2
(
7
), pp.
457
60
.10.1038/nmat924
5.
Rothstein
,
J. P.
,
2010
, “
Slip on Superhydrophobic Surfaces
,”
Annu. Rev. Fluid Mech.
,
42
(
1
), pp.
89
109
.10.1146/annurev-fluid-121108-145558
6.
Yahui
,
X.
,
Pengyu
,
L.
,
Hao
,
L.
, and
Huiling
,
D.
,
2016
, “
Underwater Superhydrophobicity: Stability, Design and Regulation, and Applications
,”
ASME Appl. Mech. Rev.
,
68
(
3
), p.
030803
.10.1115/1.4033706
7.
Rahn
,
H.
, and
Paganelli
,
C. V.
,
1968
, “
Gas Exchange in Gas Gills of Diving Insects
,”
Respiration Physiol.
,
5
(
1
), pp.
145
164
.10.1016/0034-5687(68)90083-2
8.
Lv
,
P.
,
Xue
,
Y.
,
Shi
,
Y.
,
Lin
,
H.
, and
Duan
,
H.
,
2014
, “
Metastable States and Wetting Transition of Submerged Superhydrophobic Structures
,”
Phys. Rev. Lett.
,
112
(
19
), p.
196101
.
9.
Lee
,
C.
, and
Kim
,
C.-J.
,
2011
, “
Underwater Restoration and Retention of Gases on Superhydrophobic Surfaces for Drag Reduction
,”
Phys. Rev. Lett.
,
106
(
1
).
10.
Dilip
,
D.
,
Jha
,
N. K.
,
Govardhan
,
R. N.
, and
Bobji
,
M. S.
,
2014
, “
Controlling Air Solubility to Maintain “Cassie” State for Sustained Drag Reduction
,”
Colloids Surf. A: Physicochem. Eng. Aspects
,
459
, pp.
217
224
.10.1016/j.colsurfa.2014.07.006
11.
Emami
,
B.
,
Hemeda
,
A. A.
,
Amrei
,
M. M.
,
Luzar
,
A.
,
,
M.
, and
Vahedi Tafreshi
,
H.
,
2013
, “
Predicting Longevity of Submerged Superhydrophobic Surfaces With Parallel Grooves
,”
Phys. Fluids
,
25
(
6
), p.
062108
.10.1063/1.4811830
12.
Lv
,
P.
,
Xue
,
Y.
,
Liu
,
H.
,
Shi
,
Y.
,
Xi
,
P.
,
Lin
,
H.
, and
Duan
,
H.
,
2015
, “
Symmetric and Asymmetric Meniscus Collapse in Wetting Transition on Submerged Structured Surfaces
,”
Langmuir
,
31
(
4
), pp.
1248
1254
.10.1021/la503465q
13.
Hensel
,
R.
,
Finn
,
A.
,
Helbig
,
R.
,
Killge
,
S.
,
Braun
,
H.-G.
, and
Werner
,
C.
,
2014
, “
In Situ Experiments to Reveal the Role of Surface Feature Sidewalls in the Cassie–Wenzel Transition
,”
Langmuir
,
30
(
50
), pp.
15162
15170
.10.1021/la503601u
14.
Søgaard
,
E.
,
Andersen
,
N. K.
,
Smistrup
,
K.
,
Larsen
,
S. T.
,
Sun
,
L.
, and
Taboryski
,
R.
,
2014
, “
Study of Transitions Between Wetting States on Microcavity Arrays by Optical Transmission Microscopy
,”
Langmuir
,
30
(
43
), pp.
12960
12968
.10.1021/la502855g
15.
Dufour
,
R.
,
,
N.
,
Carlier
,
J.
,
Campistron
,
P.
,
Nassar
,
G.
,
Toubal
,
M.
,
Boukherroub
,
R.
,
Senez
,
V.
,
Nongaillard
,
B.
, and
Thomy
,
V.
,
2013
, “
Acoustic Tracking of Cassie to Wenzel Wetting Transitions
,”
Langmuir
,
29
(
43
), pp.
13129
13134
.10.1021/la402481b
16.
,
J.
,
Karamanis
,
G.
,
Kirk
,
T.
, and
Hodes
,
M.
,
2017
, “
One-Dimensional Analysis of Gas Diffusion-Induced Cassie to Wenzel State Transition
,”
ASME J. Heat Transfer-Trans. ASME
,
139
(
12
), p.
122006
.10.1115/1.4036600
17.
Goswami
,
S.
,
Klaus
,
S.
, and
Benziger
,
J.
,
2008
, “
Wetting and Absorption of Water Drops on Nafion Films
,”
Langmuir
,
24
(
16
), pp.
8627
8633
.http://www.agcce.com/cytop-technical-information/
18.
Epstein
,
P. S.
, and
Plesset
,
M. S.
,
1950
, “
On the Stability of Gas Bubbles in Liquid-Gas Solutions
,”
J. Chem. Phys.
,
18
(
11
), pp.
1505
1509
.10.1063/1.1747520
19.
Oberkampf
,
W. L.
,
1976
, “
Domain Mappings for the Numerical Solution of Partial Differential Equations
,”
Int. J. Numer. Methods Eng.
,
10
(
1
), pp.
211
223
.10.1002/nme.1620100116
20.
Furzeland
,
R. M.
,
1977
, “
The Numerical Solution of Two-Dimensional Moving Boundary Problems Using Curvilinear co-Ordinate Transformations
,” Brunel University Mathematics Technical Papers Collection.
21.
John
,
C.
,
1987
,
Free and Moving Boundary Problems
,
Clarendon Press
,
New York/Oxford, UK
22.
Trefethen
,
L. N.
,
2000
,
Spectral Methods in MATLAB
, Vol.
10
,
SIAM, Society for Industrial and Applied Mathematics
,
.
23.
Weidong
,
G.
,
Gérard
,
L.
, and
Ranga
,
N.
,
2012
,
The Application of the Chebyshev-Spectral Method in Transport Phenomena
,
68
,
Springer
,
New York; Berlin
.
24.
Canuto
,
C.
,
Hussaini
,
M. Y.
,
Quarteroni
,
A.
, and
Zang
,
T. A.
,
1988
,
Spectral Methods in Fluid Dynamics
,
Springler-Verlag
,
New-York
.
25.
Sander
,
R.
,
2016
,
Henry's Law Constants
,
National Institute of Standards and Technology
,
Gaithersburg, MD
, accessed June 29, 2016, http://webbook.nist.gov
26.
Ferrell
,
R. T.
, and
Himmelblau
,
D. M.
,
1967
, “
Diffusion Coefficients of Nitrogen and Oxygen in Water
,”
J. Chem. Eng. Data
,
12
(
1
), pp.
111
115
.10.1021/je60032a036
27.
Lemmon
,
E. W.
,
McLinden
,
M. O.
, and
Friend
,
D. G.
,
Thermophysical Properties of Fluid Systems
,
National Institute of Standards and Technology
,
Gaithersburg, MD
, accessed June 29, 2016, http://webbook.nist.gov
28.
Game
,
S.
,
Hodes
,
M.
,
Keaveny
,
E.
, and
Papageorgiou
,
D.
,
2017
, “
Physical Mechanisms Relevant to Flow Resistance in Textured Microchannels
,”
Phys. Rev. Fluids
,
2
(
9
), p.
094102
.10.1103/PhysRevFluids.2.094102