1. Introduction
The behaviour of fluids when the dominant force is surface tension underpins many fundamental physical and industrially valuable processes, including microfluidics and inkjet printing (Yang, Yang & Hong Reference Yang, Yang and Hong2005; He et al. Reference He, Yang, Qin, Wen and Zhang2017; Calver et al. Reference Calver, Gaffney, Walsh, Durham and Oliver2020); directional transport of liquids in biological processes (Prakash, Quéré & Bush Reference Prakash, Quéré and Bush2008; Zheng et al. Reference Zheng, Bai, Huang, Tian, Nie, Zhao, Zhai and Jiang2010; Ju et al. Reference Ju, Bai, Zheng, Zhao, Fang and Jiang2012; Comanns et al. Reference Comanns, Buchberger, Buchsbaum, Baumgartner, Kogler, Bauer and Baumgartner2015; Xu & Jensen Reference Xu and Jensen2017; Bhushan Reference Bhushan2019); engineering applications such as water harvesting (Brown & Bhushan Reference Brown and Bhushan2016; Xu et al. Reference Xu, Lin, Zhang, Shi and Zheng2016; Li et al. Reference Li, Zhou, Li, Che, Yao, McHale, Chaudhury and Wang2017); and the behaviour of fluids in low-gravity situations (Passerone Reference Passerone2011). Moreover, from a purely theoretical point of view, such systems have been known to exhibit a plethora of complex behaviours associated with contact-line dynamics, such as contact-angle hysteresis (Dussan Reference Dussan1979; Gao & McCarthy Reference Gao and McCarthy2006; Eral, ’t Mannetje & Oh Reference Eral, ’t Mannetje and Oh2013) and the ‘stick-slip’ phenomenon (Picknett & Bexon Reference Picknett and Bexon1977; Bourges-Monnier & Shanahan Reference Bourges-Monnier and Shanahan1995; Shanahan Reference Shanahan1995; Stauber et al. Reference Stauber, Wilson, Duffy and Sefiane2014). However, the behaviour that we study here is that of a very simple static system, which forms the ‘basic state’ for many of these dynamical problems. Specifically, we seek to quantify and describe the sensitivity of menisci in confined channels to imperfections in geometry. Understanding such sensitivity is important in the industrial and biological processes described above, where natural or manufactured surfaces are in general not perfectly smooth. Indeed, the sensitivity of microfluidic devices to small imperfections has hampered their usefulness in an industrial setting (Zhou et al. Reference Zhou, Khodakov, Ellis and Voelcker2012; Calver et al. Reference Calver, Gaffney, Walsh, Durham and Oliver2020). On the other hand, the structuring of channels by parallel grooves has also been shown to improve the efficacy of microfluidic devices: for example, railed microfluidic channels have been used to create superhydrophobic surfaces (Yoshimitsu et al. Reference Yoshimitsu, Nakajima, Watanabe and Hashimoto2002; Emami et al. Reference Emami, Hemeda, Amrei, Luzar, Gad-el Hak and Vahedi Tafreshi2013), to guide and assemble microstructures inside fluidic channels (Chung et al. Reference Chung, Park, Shin, Lee and Kwon2008), and in primary cell culture technology to control the deposition and location of cells within microfluidic devices (Khademhosseini et al. Reference Khademhosseini, Yeh, Jon, Eng, Suh, Burdick and Langer2004; Howell et al. Reference Howell, Mott, Fertig, Kaplan, Golden, Oran and Ligler2005; Khademhosseini et al. Reference Khademhosseini, Yeh, Eng, Karp, Kaji, Borenstein, Farokhzad and Langer2005; Park et al. Reference Park, Toner, Yarmush and Tilles2006; Lee, Hung & Lee Reference Lee, Hung and Lee2007; Manbachi et al. Reference Manbachi, Shrivastava, Cioffi, Chung, Moretti, Demirci, Yliperttula and Khademhosseini2008; Khabiry et al. Reference Khabiry, Chung, Hancock, Soundararajan, Du, Cropek, Lee and Khademhosseini2009; Mobasseri et al. Reference Mobasseri, Faroni, Minogue, Downes, Terenghi and Reid2015). Stone, Stroock & Ajdari (Reference Stone, Stroock and Ajdari2004) also provide a general discussion of the role of channel geometry in controlling fluids in microchannels. Anna (Reference Anna2016) and Ajaev & Homsy (Reference Ajaev and Homsy2006) give a more general discussion of the modelling and applications of drops and bubbles in microchannels and confined channels.
It has been known since Wenzel (Reference Wenzel1936) that menisci in channels are sensitive to imperfections in the channel geometry. Wenzel's work demonstrating the effect of wall roughness on contact-line wettability using simple energy conservation arguments was further extended for porous surfaces by Cassie & Baxter (Reference Cassie and Baxter1944) and Cassie (Reference Cassie1948). Quasi-static effects such as contact-line hysteresis (at finite microscopic contact angle) and the stick-slip phenomenon were first observed by Johnson & Dettre (Reference Johnson and Dettre1964), who studied the wettability of a drop on a rough surface where the roughness was in the form of concentric circular grooves. By moving the contact line very slowly over the obstacles, they observed multiple axisymmetric equilibria. Huh & Mason (Reference Huh and Mason1977) studied surfaces with more complex roughness, including cross, hexagonal and radial grooves. They used the linearised Young–Laplace equation to find a relationship between the contact angle and hysteresis/stick-slip behaviour. Surfaces with periodic roughness (Cox Reference Cox1983) and random roughness (Jansons Reference Jansons1985) were also found to induce contact-angle hysteresis and stick-slip behaviour in the limit of zero capillary number. More recently, rough surfaces with Gaussian-type defects have been shown to influence significantly contact-line dynamics in the contexts of droplet spreading (Espín & Kumar Reference Espín and Kumar2015), droplet sliding (Park & Kumar Reference Park and Kumar2017), and droplet evaporation and imbibition (Pham & Kumar Reference Pham and Kumar2017, Reference Pham and Kumar2019). Jansons (Reference Jansons1985) made the further observation that the location of the contact line influenced its future position, leading to irreversibility of the wetting process. Stick-slip behaviour can also be induced by defects caused by changes in wettability or temperature (Ajaev, Gatapova & Kabov Reference Ajaev, Gatapova and Kabov2020).
Concus & Finn (Reference Concus and Finn1969), Fowkes & Hood (Reference Fowkes and Hood1998) and Reyssat (Reference Reyssat2014) showed that for static liquid–vapour interfaces in a wedge, the existence of solutions depends on the angle of the wedge and the contact angle between the liquid–vapour and solid–liquid interfaces. Instead of imposing perturbations on the channel walls geometrically, it is also possible to perturb the meniscus by changing the contact angle locally on the upper and lower channel walls. Boruvka & Neumann (Reference Boruvka and Neumann1978) considered a meniscus in contact with a stripwise heterogeneous wall in which each strip has a different equilibrium contact angle. They showed that locally near the wall, the contact-line curvature becomes unbounded at the point where the contact angle changes. The jump in contact angle is analogous to a ridge or groove perturbation with corners; here the displacement of the contact line can be unbounded in some circumstances (Concus & Finn Reference Concus and Finn1969; Weislogel & Lichter Reference Weislogel and Lichter1998; King, Ockendon & Ockendon Reference King, Ockendon and Ockendon1999). In what follows, we consider smooth (differentiable everywhere) wall perturbations with no sharp corners.
At low flow speeds, the motion of drops and bubbles in confined devices follows a series of near-equilibrium configurations and thus the static problem discussed here can also be used to provide insight into the effect of wall roughness on these problems. Channel imperfections are known to have a significant effect on the set of observable stable solutions for air fingers and bubbles propagating in a Hele-Shaw channel. This effect has been seen with a cusp at the tip of a bubble/finger created by a needle (Hong & Family Reference Hong and Family1988); a tiny bubble at the tip of a bubble/finger (Maxworthy Reference Maxworthy1986); anisotropy by etching of the plates (Ben-Jacob et al. Reference Ben-Jacob, Godbey, Goldenfeld, Koplik, Levine, Mueller and Sander1985; Dorsey & Martin Reference Dorsey and Martin1987); and channel occlusions (Hazel et al. Reference Hazel, Pailha, Cox and Juel2013; Thompson, Juel & Hazel Reference Thompson, Juel and Hazel2014). Wall roughness has also been seen to affect the ‘tip-splitting’ effect seen by propagating interfaces (Tabeling, Zocchi & Libchaber Reference Tabeling, Zocchi and Libchaber1987; Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2016, Reference Franco-Gómez, Thompson, Hazel and Juel2018). It is unclear so far how the wall roughness contributes to the selection of a particular set of solutions.
In this study, we consider a static liquid–vapour interface in a rectangular channel and introduce imperfections to the upper and lower walls in the form of narrow ridges and grooves running the length of the channel. We are interested in how the meniscus shape changes due to the perturbations and, in particular, how the perturbations displace the contact line (defined as the intersection of the meniscus with the channel walls). We consider two classes of perturbations: those that change the volume of the channel, and those that preserve it. A key result is that small-amplitude perturbations that change the volume of the channel induce a change in the pressure difference across the meniscus, and thus change the mean curvature of the meniscus. This change is a long-range effect that is felt along the whole contact line.
In § 2 we present the governing nonlinear Young–Laplace equation and boundary conditions. In § 2.1 we derive and solve a linear model to find the shape of the meniscus for perturbations of small amplitude. The linearised problem shows that the change in mean curvature of the meniscus (i.e. the change in pressure difference over the meniscus) due to small perturbations is given by the Helmholtz equation. In § 3, we derive an explicit expression for the change in pressure difference that is directly proportional to the integral of the perturbations around the contact line, which corresponds to the change in channel volume. We also solve for the fully nonlinear mean curvature of the meniscus numerically using Surface Evolver (Brakke Reference Brakke1992).
We present results for both the linear and nonlinear models in § 4. We show that for channel-volume-changing perturbations, the meniscus away from the local perturbation matches onto a catenoid that must have the same constant mean curvature as the meniscus, which can be worked out a priori from the boundary data. We also show that perturbations that are offset from each other on the two walls – for example, forming weakly corrugated channels – can be used to engineer patterns in the contact line because each perturbation causes a deflection of both the upper and lower contact lines. The deflection mechanism acts differently depending on whether the total perturbation is channel-volume-changing or channel-volume-preserving. We conclude with a short discussion in § 5.
2. Model
Consider a static liquid–vapour interface having uniform mean curvature in a rectangular channel with non-dimensional edge lengths 
$2L$, 
$2W$ and 
$1$ in the 
$x$, 
$y$ and 
$z$ directions, respectively; see figure 1. We assume that the channel is sufficiently small that gravitational effects can be ignored, i.e. that the height of the horizontal channel is much smaller than the capillary length scale 
$\sqrt {\gamma _{LV}/\rho g}$, where 
$\gamma _{LV}$ is the liquid–vapour surface tension, 
$\rho$ is the liquid density, and 
$g$ is the acceleration due to gravity. The contact angle between the liquid–vapour and solid–liquid interfaces on the upper and lower walls of the channel is fixed to be 
$\phi$, where 
$0 \leq \phi < {\rm \pi}/2$. We impose that the meniscus meets the side walls at 
$y = \pm W$ normally, with contact angle 
${\rm \pi} /2$, so that in the absence of wall perturbations, the interface takes the shape of an arc of a cylinder. To describe this, we adopt cylindrical polar coordinates 
$(r, \theta, y)$, with 
$r=0$ fixed at the centre of curvature of the unperturbed meniscus. We define the maximum value of 
$\theta$, which specifies the contact-line location in the unperturbed state, by 
${\tilde {\theta }} = {\rm \pi}/2 - \phi$ (see figure 1). Then the base state is given by 
$r= R \equiv 1/(2\sin {\tilde {\theta }})$, for 
$\theta \in [-{\tilde {\theta }}, {\tilde {\theta }}]$ and 
$y \in [-W, W]$. Cartesian and polar coordinates are related by 
$(x, y, z) = (x_0 + r \cos \theta, y, r\sin \theta )$, where 
$x_0$ is defined by the volume of liquid 
$V_L$ in the channel, as
\begin{equation} x_0 = 2L - \frac{R}{2}\cos{\tilde{\theta}} - \frac{V_L}{2W} -R^2 {\tilde{\theta}}. \end{equation}
We denote the dimensionless pressure of the liquid phase (scaled on surface tension over channel depth) to be 
$p_L$, relative to zero pressure in the gas phase, so that in the unperturbed configuration 
$p_L=-1/R$.
A static liquid–vapour meniscus in a rectangular channel 
$0 \leq x \leq 2L$, 
$-W \leq y \leq W$ and 
$-\tfrac {1}{2} + B_-(y) \leq z \leq \tfrac {1}{2} + B_+(y)$. The shape of the meniscus is described using a cylindrical polar coordinate system 
$(r, \theta, y)$ centred at 
$(x, y, z) = (x_0, 0, 0)$ where 
$x_0$ is fixed by the volume of liquid in the channel. The meniscus is located in the channel such that the unperturbed state has 
$r=R=1/(2\cos \phi )$, with 
$\phi$ the solid–liquid contact angle on the upper and lower walls. The solid–liquid contact angle on the side walls is 
${\rm \pi} /2$.
The upper (
$+$) and lower (
$-$) walls of the channel are then perturbed so that they are described by 
$z = \pm \tfrac 1 2 + B_{\pm }(y)$. The perturbations take the form of ridges and grooves that could be created, for example, using a pulsed laser (Xing et al. Reference Xing, Deng, Lian, Zhang, Zhang and Zhao2014), or using moulded fabrication with standard soft lithography techniques (Chung et al. Reference Chung, Park, Shin, Lee and Kwon2008). We do not assume contact-angle hysteresis, but we assume that surfaces are homogeneous and smooth, allowing us to address interactions between the wall perturbations and the meniscus at the microscopic level. We non-dimensionalise all surface tensions on the liquid–vapour surface tension 
$\gamma _{LV}$ so that the meniscus has unit surface tension.
We specify the interface location relative to the base state so that the surface of the meniscus is described in terms of the unknown radial perturbation 
$F(y, \theta )$ and the angular change in location of the contact lines on the upper and lower walls 
$\varPhi _\pm (y)$:
\begin{equation} r = R + F(y, \theta), \quad \theta \in [ -{\tilde{\theta}} + \varPhi_-(y), {\tilde{\theta}} + \varPhi_+(y)], \quad y \in [{-}W, W]. \end{equation}
We measure the pressure difference across the meniscus as the liquid pressure minus the gas pressure, where without loss of generality we assume that the gas pressure is zero. Thus we write 
$\Delta p = -R^{-1} + (p_L)_p$, where 
$(p_L)_p$ is the change in pressure of the liquid phase due to the channel perturbations. This is constrained by the requirement that the volume of liquid 
$V_L$ does not change under any perturbation to the channel geometry. Then, defining the unit normal 
$\boldsymbol {n}$ of the meniscus to point into the liquid phase, as shown in figure 1, the equilibrium state is specified by the Young–Laplace equation, relating the uniform mean curvature of the interface to the pressure difference across the meniscus 
$\Delta p$ as
\begin{equation} \Delta p ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{n}\vert_{r=R+F}= {-\frac{1}{\varLambda} + \left(\frac{(R+F)F_y}{\varLambda}\right)_y+ \frac{1}{R+F}\left(\frac{F_\theta}{\varLambda}\right)_\theta,} \end{equation}
where 
$\varLambda \equiv \sqrt {(R+F)^2(1+F_y^2)+F_\theta ^2}$.
We solve the Young–Laplace equation (2.3) subject to the volume constraint and the following boundary conditions. First, we constrain the contact line to lie on the perturbed channel walls, so that
\begin{equation} z = \left(R + F(y, \theta) \right) \sin\theta ={\pm} \tfrac 1 2 + B_{{\pm}}(y) \quad \mbox{at } \theta ={\pm} {\tilde{\theta}} + \varPhi_{{\pm}}(y). \end{equation}
Second, we impose a fixed contact angle 
$\phi$ through the geometrical argument that if 
$\boldsymbol {v}_{\pm }$ are unit normals to the upper and lower channel walls pointing out of the channel (as shown in figure 1), then
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{v}_{{\pm}} = \cos\phi \quad \mbox{at} \ \theta ={\pm} {\tilde{\theta}} + \varPhi_{{\pm}}(y), \quad \mbox{where } \boldsymbol{v}_{{\pm}} ={\pm} \frac{- B_{{\pm}}'(y)\,\hat{\boldsymbol{y}}+ \hat{\boldsymbol{y}}}{\sqrt{1 + B_{{\pm}}'(y)^2}}. \end{equation}Finally, the meniscus meets the side walls normally, so that
\begin{equation} F_y({\pm} W, \theta) = 0. \end{equation}2.1. Linear model
We linearise the problem when wall perturbations are small relative to the channel depth by writing 
$B_{\pm }(y) = \epsilon \,b_{\pm }(y)$, where 
$\epsilon$ is defined as the maximum amplitude of the perturbation, and 
$b_{\pm }(y) = O(1)$ as 
$\epsilon \to 0$. We use the parametrisation of the interface given in the nonlinear model (2.2), with the assumption that the perturbation to the radius and the change in location of the contact lines are also 
$O(\epsilon )$, so that 
$F(y, \theta ) = \epsilon \,f(y, \theta )$ and 
$\varPhi _\pm (y) = \epsilon \,\varTheta _\pm (y)$. Then the interface is parametrised as
\begin{equation} r = R + \epsilon\,f(y, \theta), \quad \theta \in \left[- {\tilde{\theta}}+ \epsilon\, \varTheta_-(y), {\tilde{\theta}} +\epsilon\, \varTheta_+(y) \right], \quad y \in [{-}W, W]. \end{equation}
We assume that the change in the pressure difference due to the perturbation is also 
$O(\epsilon )$, so that 
$(p_L)_p = \epsilon p$ and 
$\Delta p = -R^{-1} + \epsilon p$.
Assuming that 
$f(y, \theta )$ and all its derivatives are 
$O(1)$ as 
$\epsilon \to 0$, after linearisation, the leading-order approximation to the Young–Laplace equation (2.3) is the Helmholtz equation
\begin{equation} \frac{1}{R^2}\,f +\frac{1}{R^2}\,f_{\theta\theta} + f_{yy} =p, \quad y \in [{-}W, W], \ \theta \in [-{\tilde{\theta}}, {\tilde{\theta}}]. \end{equation}The boundary condition (2.4) constraining the contact line to lie on the upper and lower channel walls becomes
\begin{equation} {f(y, \pm {\tilde{\theta}})\sin{\tilde{\theta}} \pm f_\theta(y,\pm {\tilde{\theta}}) \cos{\tilde{\theta}} ={\pm} b_\pm(y)}, \end{equation}whilst the leading-order approximation to boundary condition (2.5) is
\begin{equation} {f_\theta(y, \pm {\tilde{\theta}}) - R \, \varTheta_\pm (y) = 0.} \end{equation}
In deriving (2.10), the neglected 
$O(\epsilon ^2)$ terms include those involving the derivative of the boundary data, 
$\epsilon ^2\,b_{\pm }'(y)$. These terms will formally become 
$O(\epsilon )$ if 
$b_{\pm }'(y) = O(\epsilon ^{-1})$; this puts an additional constraint on the boundary data for the linearised model to be valid. In particular, this constraint suggests that we cannot use a linearised model for very sharp perturbations with large gradients regardless of their amplitude; this will be discussed in § 3.2. The final boundary condition (2.6) becomes
\begin{equation} f_y({\pm} W, \theta) = 0. \end{equation}
The problem is closed with the condition that the volume of liquid 
$V_L$ is invariant with respect to changes in channel volume. This leads to the condition
\begin{equation} \int_{{-}W}^{W} \int_{-{\tilde{\theta}}}^{{\tilde{\theta}}} f(y, \theta) \, \mathrm{d}\theta \,{{\rm d} y} = \left(\frac{V_L}{2RW} + \frac{{\tilde{\theta}}}{2\sin{\tilde{\theta}}}-\frac{\cos{\tilde{\theta}}}{2}\right) \int_{{-}W}^{W} \left( b_+(y) -b_-(y)\right) {{\rm d} y}. \end{equation}
Finally, the linearised contact-line displacement on the upper and lower walls 
$x_\pm (y)$ is found by Taylor expanding the solution for 
$x = r \cos \theta$ at 
$\theta = \pm {\tilde {\theta }} + \epsilon \,\varTheta _\pm (y)$:
\begin{align} x_\pm(y) = x_0 + R \cos{\tilde{\theta}} + \epsilon\,x_p^\pm(y) + O(\epsilon^2), \quad x_p^\pm(y) = f(y, \pm {\tilde{\theta}})\cos {\tilde{\theta}} \mp f_\theta(y, \pm {\tilde{\theta}}) \sin{\tilde{\theta}}. \end{align} For the specific case of zero contact angle (
$\phi =0$), when the meniscus meets the walls tangentially, there is no contribution to boundary condition (2.10) at 
$O(\epsilon )$. An expansion to powers of 
$O(\epsilon ^2)$ is needed to obtain (2.10), as explained in Appendix A.
3. Methods
3.1. Pressure-volume relationship
We show that small-amplitude perturbations that change the volume of the channel must cause the mean curvature of the meniscus to change. Noting the self-adjointness of the Helmholtz operator and boundary conditions (2.8)–(2.11), we derive in Appendix B an explicit expression for this 
$O(\epsilon )$ pressure difference at the meniscus, and find that it has a similar dependence on the boundary data as the induced volume change,
\begin{equation} p = \frac{1}{4 W R^2 \sin ({\tilde{\theta}})} \int_{{-}W}^{W} \left( b_+(y) - b_-(y)\right) {{\rm d} y}. \end{equation}
Thus in the linear problem with channel-volume-changing perturbations (for which 
$\int _{-W}^{W} ( b_+(y) - b_-(y)) \, {{\rm d} y} \neq 0$), the forcing for the Helmholtz equation (2.8) can be deduced a priori from the volume change encoded in the boundary data.
3.2. Solutions for Gaussian perturbations
We now restrict our attention to Gaussian perturbations of the form
\begin{equation} {B_\pm(y) = a^\pm }\epsilon \exp({-(y-y_c^\pm)^2/s^2}), \quad \mbox{so}\quad b_\pm(y) = a^\pm \exp({-(y-y_c^\pm)^2/s^2}), \end{equation}
where 
$y_c^\pm$ and 
$s$ control the location and width of the perturbation, respectively, and the prefactor 
$a^\pm$, which can take value 
$+1$ or 
$-1$ on either wall, determines the orientation of the perturbation, i.e. whether it is a ridge or a groove. For the purposes of illustration, we assume that the perturbation on the lower wall is a ridge so that 
$a^-=1$. Then after fixing 
$s$, we consider two specific types of geometry: channel-volume-preserving configurations with 
$\int _{-W}^{W} ( b_+(y) - b_-(y)) \, {{\rm d} y} = 0$ (corresponding to a groove on the upper wall); and channel-volume-changing configurations with 
$\int _{-W}^{W} ( b_+(y) - b_-(y)) \, {{\rm d} y} < 0$ (corresponding to a ridge on the upper wall). The former preserve the pressure of the liquid phase, whereas the latter, which decrease the volume of the channel, cause an increase in the magnitude of the liquid pressure. If 
$y_c^\pm =0$, then the channel-volume-preserving and channel-volume-changing configurations are mirror-antisymmetric and mirror-symmetric, respectively, about 
$z=0$.
We solve the full nonlinear problem to find the shape of the interface using the open-source software Surface Evolver (Brakke Reference Brakke1992), which uses a gradient-descent method to converge to a surface with minimum energy from a given initial guess and subject to constraints to enforce the boundary conditions and the volume condition; for details, see Appendix C. Meanwhile, we solve the linear problem (2.8)–(2.12) using second-order-accurate finite differences; see Appendix D for more details. As seen in § 2.1, the linear model can be expected to break down when 
$b_\pm '(y) = O(\epsilon ^{-1})$. For the Gaussian boundary data (3.2), this occurs if 
$s^2 =O(\epsilon )$, which limits how narrow the perturbation can be.
3.3. Analytic solution for symmetric perturbations
For the special case of aligned perturbations (
$y_c^\pm =0$), we can obtain an analytic solution to the linear problem (2.8)–(2.12) via separation of variables, with a Fourier discretisation across the width of the channel in the 
$y$ direction because of the 
$y$ dependence of the boundary conditions on the upper and lower walls. Denoting mirror antisymmetric (channel-volume-preserving) and mirror-symmetric (channel-volume-changing) solutions by ‘MAS’ and ‘MS’ subscripts, respectively, the series solution is given by
\begin{equation} f_{{MAS}\atop{MS}}(y, \theta) = R^2 p+ \sum_{n=0}^\infty A_n^{{MAS}\atop{MS}} \left( \mathrm{e}^{\lambda_n\theta } \mp \mathrm{e}^{-\lambda_n\theta } \right) \cos\left(\frac {n{\rm \pi} y }{W}\right), \end{equation}
where the exponential coefficient 
$\lambda _n$ is given by
\begin{equation} \lambda_n = \sqrt{\left(\frac{n{\rm \pi} R}{W}\right)^2 -1}. \end{equation}
Thus there is a critical value of 
$n$, specific to the contact angle (through 
$2R \cos \phi =1$) and the width of the channel, at which the sum switches from having oscillatory dependence in 
$\theta$ to exponential dependence in 
$\theta$. The coefficients 
$A_n^{{MAS}\atop {MS}}$ are given by
$$\begin{gather} A_n^{{MAS}\atop{MS}} = \frac {\pm a_n}{ \left(\sin {\tilde{\theta}} + \lambda_n\cos {\tilde{\theta}} \right) {{\rm e}^{\lambda_n{\tilde{\theta}}}} \mp \left( \sin {\tilde{\theta}} -\lambda_n\cos {\tilde{\theta}}\right) {{\rm e}^{-\lambda_n{\tilde{\theta}}}}}\quad \begin{cases} n \geqslant 0 \ (\text{MAS}), \\ n \geqslant 1 \ (\text{MS}), \end{cases} \end{gather}$$
$$\begin{gather}a_0^{{MAS}} = \frac {1}{2W} \int_0^{W} b_-(y) \, {{\rm d} y}, \quad a_n = \frac {1}{W} \int_0^{W} b_-(y) \cos\left(\frac {n{\rm \pi} y }{W}\right) {{\rm d} y}, \quad n \geqslant 1. \end{gather}$$
The coefficient 
$A_0^{{MS}}$ is found using the volume condition (2.12) to be
\begin{equation} A_0^{{MS}} = \frac{\left(\dfrac{V_L}{2RW} + \dfrac{{\tilde{\theta}}}{2\sin{\tilde{\theta}}} - \dfrac{\cos{\tilde{\theta}}}{2}\right) \displaystyle \int_{{-}W}^{W} \left( b_+(y) - b_-(y)\right) {{\rm d} y} - 4 W R^2 {\tilde{\theta}} \,p}{4 W \sin{\tilde{\theta}}}. \end{equation}For the Gaussian boundary data (3.2), the coefficients of the convergent series are defined by
\begin{equation} a_n = \frac {s\sqrt{ {\rm \pi}}} {2W} \exp\left(-\frac {s^2} 4 \left(\frac {n{\rm \pi} }{W} \right)^2\right)\mbox{Re} \left\{ \mathrm{i}\,\mathrm{erfi} \left( -\frac{\mathrm{i} W}{s} - \frac {n {\rm \pi}s}{2W} \right)\right\}, \quad n \geqslant 1. \end{equation}
The function 
$\mathrm {erfi}(z) = -\mathrm {i}\,\mathrm {erf}(\mathrm {i} z)$ is the imaginary error function so that for real 
$u$ and 
$v$, 
$\mathrm {i}\,\mathrm {erfi}(-\mathrm {i} u+v) = \mathrm {erf}(u+\mathrm {i} v )$.
Numerical solutions below are obtained by truncating (3.3) at 
$n=n_c$ such that terms with coefficients smaller than 
$10^{-16}$ were discarded.
4. Results
We present results for the displacement of the static meniscus and the contact line induced by Gaussian perturbations (3.2). We first assume that the perturbations are aligned so that 
$y_c^\pm = 0$.
4.1. Aligned perturbations
Figure 2 shows ‘baseline’ linear solutions 
$f(y, \theta )$ for the two prototype channel-volume- preserving and -changing configurations, together with displacement of the upper and lower contact lines due to the perturbation, 
$x_p^\pm (y)$ (as given in (2.13a,b)), computed using the series solution (3.3). In the channel-volume-preserving case (figure 2a), the response of the meniscus and contact line is localised around the perturbations in the 
$y$ direction, whereas in the channel-volume-changing case (figure 2b), the perturbations induce a larger-amplitude response that is felt across the entire depth and width of the channel. In the former case, the contact-line shape appears to mirror the curvature of the wall perturbation, but it is smoother in the latter case. Thus a small ridge or groove placed in the centre of the channel can cause non-local bending of the contact line through its impact on the pressure field (3.1). Note that we can build new small-amplitude solutions as linear combinations of the two solutions shown, and can therefore describe the behaviour of the contact line as the perturbations are varied between the two configurations.
The displacement of the upper and lower contact lines due to the perturbation 
$x_p^\pm (y)$ (black solid lines, left axis) and change in meniscus shape 
$f(y, \theta )$ (heat map) due to Gaussian perturbations 
$b_\pm (y)$ proportional to 
$b(y) = \exp (-y^2/0.75)$ (red dashes, right axis), for a channel with half-width 
$W = 5$ and contact angle 
$\phi = 15^{\circ }$. As in figure 1, the liquid and vapour sides of the contact-line displacement are shown with blue and white shading, respectively, for (a) channel-volume-preserving perturbation (
$b_- = b_+ = b$); and (b) channel-volume-changing perturbation (
$b_- = -b_+ = b$). The heat maps denote the change in radius 
$r$ due to the perturbations, with green denoting no change in the meniscus location. Positive values of 
$f$ show the meniscus extending into the liquid phase.
The effect of varying perturbation amplitude and width on the contact-line displacement is presented in figure 3. Both volume-preserving (figure 3a) and volume-changing (figure 3b,c) perturbations result in a local displacement of the contact line in the centre of the channel that depends on the shape of the perturbations, with narrower or larger-amplitude perturbations eliciting a greater displacement. Linear and nonlinear predictions of the upper contact-line displacement 
$x_p^+(y)$ show good agreement for perturbations that are approximately 
$1\,\%$ of the channel depth (figure 3a,b). (Numerical evidence in figure 10 below suggests that the contact-line displacement increases like 
$\log (1/s)$ for finite contact angle as the perturbation width 
$s$ becomes very small, when the linear model breaks down.) Larger perturbations, of approximately 
$10\,\%$ of channel depth, are shown in figure 3(c) for the volume-changing case. Here a greater discrepancy between linear and nonlinear predictions is evident, although the shape of the contact-line perturbations at the edges of the channel is accurately described by the linear model. A long-wave analysis for 
$s \ll W$ shows that the solution in the far-field has a quadratic dependence on 
$y$ in the form 
$f_A(y, \theta ) = C_1(\theta )\,(W-|y|)^2+C_2(\theta )$. We substitute this ansatz into the Helmholtz equation (2.8) and solve the resulting coupled ODEs for 
$C_1(\theta )$ and 
$C_2(\theta )$ with modified ‘far-field’ boundary conditions. These are obtained by setting 
$b_\pm (y) = 0$ in (2.9) since the amplitude of Gaussian perturbations 
$b_\pm (y)$ is negligibly small for 
$y \sim W$ and 
$\vert y-y_c\vert \gg s$. Then substituting the solution for 
$f_A(y, \theta )$ in (2.13a,b) shows that the approximation to the displacement of the contact line far from the perturbations is given by
\begin{equation} \hat x_p^+ (y; {\tilde{\theta}}) \approx \left({\frac {p\sin ( {\tilde{\theta}} ) }{\cos ( {\tilde{\theta}} ) \sin ( {\tilde{\theta}} ) +{\tilde{\theta}}}}\right)(W-|y|)^2 + C, \end{equation}
where the translational constant 
$C$ cannot be found using the volume condition and instead is found empirically by comparison with the value of the full solution (3.3) at 
$y= \pm W$. This linear ‘far-field solution’, which is valid when 
$\vert y-y_c\vert \gg s$, is shown in figure 3(c) and gives an excellent fit to the nonlinear data.
The displacement of the upper contact line in channels of half-width 
$W=5$, for perturbations with amplitude 
$\epsilon = 0.01$ (a,b) and 
$\epsilon =0.1$ (c), for channel-volume-preserving perturbations (a) and channel-volume-changing perturbations (b,c). The contact angle is 
$\phi = 85^{\circ }$. Panels (a,b) show wall perturbations of varying width, with darker colours indicating wider perturbations, with 
$s^2$ varying from 
$0.75$ (black), through 
$0.5$ and 
$0.25$, to 
$0.1$ (light grey); (c) shows 
$s^2 = 0.75$ only. Solid lines denote the displacement calculated via the linear solution, i.e. 
$x_p^+(y)$ from (2.13a,b), for all values of 
$s$; circles (thicker lines) denote the nonlinear displacement 
$(x_{{cl}}-x_0-R\cos {\tilde {\theta }})/\epsilon$, where 
$x_{{cl}}$ is the contact-line data computed in Surface Evolver for 
$s^2 \geqslant 0.25$. Inset: the nonlinear lower contact-line displacement for volume-preserving perturbations, for 
$s^2=0.75$. The red dots in (c) denote the quadratic linear far-field solution (4.1), with 
$C \approx 0.18$.
Recalling from (3.1) that 
$p=O(W^{-1})$, the linear far-field solution (4.1) for volume-changing perturbations suggests that if the channel is sufficiently wide, so that 
$W \sim O(\epsilon ^{-1})$, then the far-field contact-line displacement could become 
$O(1)$, violating the small displacement assumption of the linear model. We therefore need to revisit the far-field quadratic approximation (4.1), which should correspond to the arc of a circle having curvature equivalent to that of a large, flat ‘pancake’ catenoid confined between unperturbed plates having the same mean curvature, i.e. the same mean pressure difference 
$\Delta p$ as the meniscus. The radius 
$R_d$ of the circular arc that matches onto the contact line is found by computing a catenoid with the same pressure difference 
$\Delta p$, as evaluated in Appendix E. Figure 4 shows how the computed far-field shape of the contact line is captured well by both the linear far-field solution (4.1) and the circular arc computed from the catenoid solution.
The nonlinear contact-line displacement for a channel of half-width 
$W=10$, with contact angle 
$\phi = 45^\circ$ and perturbations to the channel wall of amplitude 
$\epsilon = 0.01$ and 
$s^2=1$. The thick black solid line denotes the nonlinear upper contact-line data computed in Surface Evolver, 
$x_{{cl}}$. The red line denotes the quadratic linear far-field solution (4.1), 
$x_0 + R \cos {\tilde {\theta }} + \epsilon \,\hat {x}_p^+(y)$, with 
$C \approx 0.18$. The blue dashes denote the arc of a circle of radius 
$R_d$, where 
$R_d$ is determined as part of the solution to the boundary-value problem described in Appendix E with 
$\Delta p \approx -1.4167$ found from Surface Evolver.
In summary, the pressure change induced by the net volume changes of the wall perturbations generates curvature of the contact line away from the perturbations, whereas other geometric features of the perturbations influence contact-line shapes locally. We therefore expect that the same far-field behaviour should exist if the perturbations are not aligned, as we shall test in the next section.
4.2. Non-aligned perturbations
We now consider perturbations that are not aligned, i.e. 
$y_c^\pm \neq 0$. We consider specifically configurations of perturbations that are sufficiently far apart to be considered as isolated perturbations.
We compute the non-aligned solutions to the linear model using second-order-accurate central finite differences (Appendix D), with step sizes 
$\Delta y = 0.05$ in the 
$y$ direction and 
$\Delta \theta = 0.03$ in the 
$\theta$ direction. Figure 5 shows solutions to the linear problem for perturbations that have been separated so that 
$b_\pm (y)$ are centred at 
$\pm y_c$. The separation of the channel-volume-preserving perturbations causes the contact line to bend away from the side wall; this deflection of the upper and lower contact lines occurs due to the presence of a perturbation on either the upper or lower wall. Isolated ridges and grooves cause the contact lines to move towards the liquid and vapour phases, respectively. In contrast, channel-volume-changing perturbations (figure 5c,d) induce non-local bending of the contact line in the far-field, which can again be described using arcs of circles.
The linear solution 
$f(y, \theta )$ and upper contact-line displacement for channel-volume-preserving (a,b) and channel-volume-changing (c,d) perturbations, for channel half-width 
$W = 5$ and contact angle 
$\phi = 85^{\circ }$. The upper and lower wall perturbations are given by 
$B_\pm (y) = 0.01 \exp ((y - y_c^\pm )^2/0.25)$. In (a,c) 
$y_c^\pm =\pm 3$, while (b,d) show contact-line displacement 
$x_p^+$ for separations varying from 
$y_c^\pm = \pm 3$ (black) through 
$\pm (1.5, 1, 0.5)$ to 
$y_c^\pm = 0$ (light grey). Solid lines denote the displacement calculated via the linear solution 
$x_p^+(y)$ from (2.13a,b); circles (thicker lines) denote the nonlinear displacement 
$(x_{{cl}}-x_0-R\cos {\tilde {\theta }})/\epsilon$, where 
$x_{{cl}}$ is the upper contact-line data computed in Surface Evolver. The pink line in (b) is the line 
$x = \alpha y$, where the slope 
$\alpha$ is given in (4.2).
We wish to understand the deflection mechanism so that we may choose perturbations to engineer specific contact-line shapes. In the channel-volume-preserving case (figure 5a,b), let 
$\alpha$ be the gradient of the contact-line displacement in the centre of the channel, between the perturbations. For perturbations of sufficiently small amplitude, 
$\alpha$ is approximately equal to the angle of deflection, i.e. the angle that the contact line makes with the horizontal. We can obtain an approximation for 
$\alpha$ by considering the solution in the neighbourhood of an isolated ridge or groove: consider the Helmholtz equation (2.8), with zero pressure difference 
$p$ (to obtain contact-line solutions with uniform gradient). Then the problem for a single isolated ridge or groove perturbation 
$b_+(y)$ on the upper wall, centred at some 
$y=y_c$, will admit a solution of the Helmholtz equation with 
$f(y, \theta ) = 0$ for 
$y_c-y \gg s$ and 
$f(y, \theta ) \approx \alpha (y-y_c) \cos \theta$ for 
$y-y_c \gg s$. Exploiting the self-adjointness of the Helmholtz operator and boundary conditions using the method given in Appendix B with a test function 
$g(\theta ) = \cos \theta$, we obtain
\begin{equation} \alpha \approx \frac{1}{R^2}\,\frac{1}{\cos {\tilde{\theta}} \sin{\tilde{\theta}} + {\tilde{\theta}}}\int_{-{\infty}}^{\infty} b_+(y) \, {{\rm d} y}, \end{equation}
where the integral is taken over the full width of the isolated perturbation 
$b_+$. Thus we anticipate that for Gaussian perturbations, the parameters that most affect the deflection will be the volume of the perturbation and the contact angle. While the linear theory allows for an 
$O(\epsilon )$ contact-line displacement in the 
$x$ direction, the displacement 
$\epsilon \alpha y$ of the deflected solution can in principle become 
$O(1)$ in a sufficiently wide channel (if 
$y-y_c = O(\epsilon ^{-1})$); thus the solution can in principle be matched to a straight meniscus for which the 
$x$ displacement is larger. Figure 6 shows the values of 
$\alpha$ found empirically, together with the theoretical prediction (4.2), for varying perturbation separation, width, amplitude and contact angle. There is excellent agreement with the theoretical predictions except for small 
$y_c^+$, i.e. as long as the perturbations are not too close together; this is expected because it violates the assumption that the perturbations can be treated as isolated ridges and grooves.
The slopes 
$\alpha$ of the upper contact-line displacement 
$x_p^+(y)$ in the centre of the channel for channel-volume-preserving perturbations. (a) Varying contact angle 
$15^\circ \leq \phi \leq 85^\circ$ with separation 
$y_c^+ = 3$, perturbation width 
$s = 0.5$, and perturbation amplitude 
$\epsilon = 0.01$. (b) Varying amplitude 
$0.001 \leq \epsilon \leq 0.065$ with 
$y_c^+ = 3$, 
$s=1$ and 
$\phi = 85$. (c) Varying separation 
$0 \leq y_c^\pm \leq 3$ with 
$s=0.5$, 
$\epsilon = 0.01$ and 
$\phi = 85$. (d) Varying square perturbation width 
$0.25 \leq s^2 \leq 1$ with 
$y_c^+ = 3$, 
$\phi = 85$ and 
$\epsilon = 0.01$. The solid lines denote the values of 
$\alpha$ calculated using (4.2). The diamonds denote the numerical values of 
$\alpha$ calculated empirically from the contact-line data.
4.3. Weakly corrugated channels
Based on the discussion above, we now consider the linear model for channels with a series of small-amplitude ridges and grooves on the upper wall to form weakly corrugated channel walls. Thus we consider perturbations of the form
\begin{equation} b_\pm(y) = \sum_{k=0}^{K} a_k^{{\pm}} \exp({-(y-y_{c_k}^\pm)^2/s^2}), \end{equation}
where 
$y_{c_k}^\pm$ are the locations of the ridges and grooves on the upper and lower walls, and 
$a_k^{\pm }=\pm 1$ depending on whether ridges or grooves are chosen. We assume that the ridges and grooves are spaced sufficiently so that we can treat each perturbation as a single isolated ridge or groove that causes deflection of the contact line in the way described above, so that we can predict the deflection angle due to each ridge and groove using (4.2). Figure 7 shows the upper and lower contact-line displacements 
$x_p^\pm (y)$ for a weakly corrugated channel with alternating ridges and grooves on each wall so that the contact lines take the shape of a letter ‘M’. The contact-line displacement for the channel-volume-preserving configuration is shown in figure 7(a); this solution describes a meniscus with zero induced mean curvature, thus the contact-line displacement is flat in the far-field and has sections of varying slope. Because each perturbation can be treated in isolation, the gradient of each slope is described by (4.2). Again, ridges push the contact line towards the liquid phase, and grooves allow the contact line to move towards the vapour phase. The contact-line slope varies smoothly, and again the shape of the contact line is affected by an obstacle on either wall so that it takes, for example, a ridge/groove on the lower wall followed by a ridge/groove on the upper wall to reverse the gradient of the contact line.
The upper and lower contact lines (left axis, solid black line), together with the perturbations (right axis, dashed red line) for (a) a channel-volume-preserving configuration and (b) a channel-volume-changing configuration in a channel of half-width 
$W=20$. The perturbations are defined by (4.3), with ridges on the lower wall at 
$y = -14, 2$, grooves on the lower wall at 
$y=-6, 10$, ridges on the upper wall at 
$y=-2, 14$ and grooves on the upper wall at 
$y=-10, 6$. The perturbation in (b) has an extra ridge on the lower wall at 
$y=0$ to make it channel-volume changing. All perturbations have width 
$s^2=0.1$ and the contact angle is 
$\phi = 85^\circ$.
In figure 7(b), we consider the same series of grooves and ridges and then add an extra ridge on the lower wall at 
$y=0$ so that the configuration is now channel-volume-changing. Qualitatively, the shape is unchanged but the change in mean curvature is now non-zero. Thus the ‘sections’ of contact line between each ridge and groove are now locally parabolic; each parabola is a local approximation to the arc of a circle that forms the contact line of a catenoid with the same mean curvature as the static meniscus.
There are, of course, a plethora of possible smoothly-varying contact-line shapes that can be made using channel-volume-changing/preserving configurations, for small-amplitude perturbations (up to approximately 10 % of the channel height). It is possible to specify these shapes a priori using just the boundary data, using either (4.2) for the required gradients in the channel-volume-preserving case, or (3.1) to deduce the pressure of the catenoid with circular contact line that matches onto the parabolas in the channel-volume-changing case, together with the known direction in which the contact line will move for either ridges or grooves.
5. Discussion
In this study, we have quantified the displacement of the contact line of a static meniscus in a rectangular channel arising from the presence of isolated ridges and grooves imposed on the channel walls. We have shown that small-amplitude perturbations that change the channel volume induce a change in the mean curvature of the meniscus, inducing long-range curvature of the contact line, via (3.1). For very wide channels, this curvature matches onto the arc of a catenoid whose radius is found by matching the pressure differences. Meanwhile, small-amplitude isolated non-aligned perturbations that do not change the channel volume generate a contact-line shape that is approximately piecewise linear. We derived an approximation to the deflection angle between adjacent linear segments (4.2), showing a dependence on the volume of the groove or ridge. This makes it possible in principle to engineer contact-line shapes by choosing the location and order of the ridges and grooves.
We validated predictions of the linearised model against fully nonlinear solutions obtained using Surface Evolver. However it remains unclear at present how the closed-form results (3.1) and (4.2), derived using the self-adjointness of the Helmholtz equation, might be extended to the nonlinear regime. While these predictions of the induced pressure and deflection angle show dependence on the volume of ridges or grooves, they mask more subtle dependence on the precise shape of the perturbations. For example, when there is no induced pressure change, the contact-line displacement near a ridge or groove mirrors approximately the curvature of the wall shape (figure 2a), which is a bounded function for the Gaussian wall perturbations chosen here. Sharper perturbations, having derivatives varying on very short length scales, can be expected to lead to dramatically different outcomes, as outlined in Appendix F. We avoided these extreme cases here by ensuring that 
$b_{\pm }(y)$ is analytic and not too narrow.
A natural extension of this study is to consider perturbations with curvature in two directions (such as isolated bumps). These too can be expected to generate long-range deflections of the contact line. However, nonlinear effects (associated with large amplitudes or sharp asperities) will likely need to be taken into account in order to capture effects such as contact-line hysteresis, arising as the contact line is moved slowly backwards and forwards over the bump. Similarly, the approach taken here could equally be extended to consider the sensitivity of the meniscus to changes in the contact angle arising from coating portions of the channel wall with suitable chemicals. In practice, however, a continuous gradient of contact angle may be much more difficult to achieve experimentally than smoothly-varying perturbations, which may appear naturally in an industrial or biological setting. The present study considered perturbing a ‘straight’ meniscus with zero Gaussian curvature; a further generalisation that merits investigation is to consider a curved base state, to accommodate contact angles at lateral walls that deviate from 
${\rm \pi} /2$.
The solution structures identified here will support future studies of gas/liquid interfaces moving at low capillary numbers through domains having isolated geometric features, be these engineered in order to achieve a specific outcome or naturally occurring roughness. We have shown that even when these features are smooth, isolated and of small amplitude, significant long-range deflections of the meniscus are possible.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Linearised problem for zero contact angle
When 
$\phi = 0$, the linearised boundary condition (2.10) vanishes, requiring expansion up to 
$O(\epsilon ^2)$. We write the interface location as
$$\begin{gather} r = R + \epsilon\,f_1(y, \theta) + \epsilon^2\,f_2(y, \theta) + O(\epsilon^3), \\ \theta \in \left[- {\tilde{\theta}}+ \epsilon\,\varTheta_{1-}(y) + \epsilon^2\,\varTheta_{2-}(y), \ {\tilde{\theta}} +\epsilon\, \varTheta_{1+}(y) + \epsilon^2\,\varTheta_{2+}(y) \right],\nonumber\end{gather}$$
where 
$R = \tfrac 1 2$. The pressure difference 
$\Delta p$ is assumed to be 
$\Delta p = -R^{-1} + \epsilon p_{1} + \epsilon ^2 p_{2}+ \cdots$. After linearising the Young–Laplace equation (2.3), the 
$O(\epsilon )$ expression gives the equation for 
$f_1$ as
\begin{equation} \frac 1 {R^2}\,f_1 + \frac 1 {R^2}\,f_{1\theta\theta} + f_{1yy} =p_{1}. \end{equation}
Similarly, the 
$O(\epsilon )$ terms in the linearisation of boundary conditions (2.4) and (2.6) give the boundary conditions on 
$f_1$ as
\begin{equation} f_1\left(y, \pm \frac {\rm \pi}2 \right) = b_{{\pm}}(y), \quad f_{1y}({\pm} W, \theta) = 0. \end{equation}
The equation relating the change in meniscus shape 
$f(y, \theta )$, to the contact-line location 
$\varTheta _{1 \pm }(y)$, can be found at 
$O(\epsilon ^2)$:
\begin{equation} f_{1\theta}\left(y, \pm \frac {\rm \pi}2\right) = R \, \varTheta_{1\pm}(y). \end{equation}
Therefore we recover exactly the conditions (2.9) and (2.10) with 
$\phi =0$. The volume constraint (2.12) and independent pressure condition (3.1) remain the same.
Appendix B. The pressure in the linear problem
For the linearised problem, we can derive an independent equation for the pressure by using the fact that the Helmholtz operator is self-adjoint. Consider the linear problem (2.8)–(2.12) and a smooth, twice-differentiable test function 
$g(\theta ): [-{\tilde {\theta }}, {\tilde {\theta }}] \to \mathbb {R}$, such that
\begin{equation} R^{{-}2} [ g'' + g] = a \in \mathbb{R}, \quad g({\pm} {\tilde{\theta}}) = \gamma_{{\pm}}, g'({\pm} {\tilde{\theta}}) = \zeta_\pm. \end{equation}
We multiply the Helmholtz equation (2.8) by 
$g$, and then integrate over the domain 
$D = [-{{\tilde {\theta }}}, {{\tilde {\theta }}}] \times [-W, W]$. Then, defining 
$\tilde {\boldsymbol {\nabla }}= (\partial _y, R^{-1}\partial _\theta )$, we obtain
\begin{equation} \int_D a f + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot} (g \tilde{\boldsymbol{\nabla}} f - f \tilde{\boldsymbol{\nabla}}g) \,\mathrm{d} A = \int_D gp \,\mathrm{d} A. \end{equation}
Rewriting the divergence terms on the left-hand side as integrals over closed curves, we then integrate and apply the boundary conditions at the side walls, (2.11), and the boundary conditions on 
$g$ to give
\begin{align} &\int_D a f \,
\mathrm{d} A + R^{{-}2} \int_{{-}W}^{W} [-
\gamma_-\,f_\theta(y, -{{\tilde{\theta}}}) + \zeta_-\,f(y,
-{{\tilde{\theta}}}) + \gamma_+\,f_\theta(y,
{{\tilde{\theta}}}) - \zeta_+\,f(y, {{\tilde{\theta}}})]
\,{{\rm d} y}\notag\\ &\quad= \int_D gp \,\mathrm{d} A.
\end{align}
We now pick a test function 
$g(\theta ) = \cos \theta$ (so that 
$a=0$) and apply the boundary conditions (2.9) on 
$f$, which leads to an independent equation for the pressure:
\begin{equation} p = \frac{1}{4 W R^2 \sin ({{\tilde{\theta}}})} \int_{{-}W}^{W} [ b_+(y) - b_-(y)] \, {{\rm d} y}. \end{equation}We also use this method to derive an approximation to the deflection angle 
$\alpha$ as given in (4.2). To derive this relationship, we again use a test function 
$g(\theta ) = \cos \theta$, but we take ‘far-field’ boundary conditions across an isolated wall perturbation of 
$f_y \rightarrow 0$ for 
$(y-y_c)/s \rightarrow -\infty$, and 
$f_y \rightarrow \alpha \cos \theta$ for 
$(y-y_c)/s \rightarrow \infty$.
Appendix C. Solving the nonlinear problem with Surface Evolver
We implement the nonlinear problem using Surface Evolver (Brakke Reference Brakke1992), which uses a gradient-descent method to iterate towards a surface of minimum energy. The energy of the triangulated surface is defined as a scalar function of all the vertex coordinates. The iterative process forces the vertices into a configuration that is closer to an energy minimum, subject to any global constraints on the surface or local constraints on the vertices. For the problem outlined in § 2, the only force acting on the liquid–vapour interface is surface tension, thus we minimise the surface energy of the meniscus. The constraints are the boundary conditions (2.4)–(2.6) together with the volume constraint.
The boundary conditions on the upper and lower channel walls, (2.5) and (2.6), are implemented by fixing the energy of the channel walls. We define the (non-dimensional) surface tension due to the presence of the solid wall as 
$\gamma _S = \gamma _{SL}/\gamma _{LV}-\gamma _{SV}/\gamma _{LV}$, where 
$\gamma _{SL}$ and 
$\gamma _{SV}$ are solid–liquid and solid–vapour surface tensions. Then if 
$\alpha _w$ is the contact angle at the solid–liquid and liquid–vapour interface, by Young's equation, 
$\gamma _S = -\cos \alpha _w$, and the wall energy is
\begin{equation} E_{{wall}} = \iint_{{wall}} -\cos\alpha_w \,\mathrm{d} S. \end{equation}
Thus to specify the boundary conditions (2.5) and (2.6), we fix the wall energies by imposing contact angles 
$\alpha _w = \phi$ on the upper and lower walls at 
$z = \pm 1/2 + B_{\pm }(y)$, and 
$\alpha _w = {\rm \pi}/2$ on the side walls at 
$y = \pm W$ (so that the energy of the side walls is zero).
In practice, for computational efficiency, only the liquid–vapour interface is triangulated and refined. Then the wall energy integral (C1) for the upper and lower walls is rewritten as a line integral using Stokes’ theorem: if 
$\boldsymbol {w} = (w_x, w_y, w_z)$ is such that 
$(\boldsymbol {\nabla } \times \boldsymbol {w})\boldsymbol {\cdot } \boldsymbol {v}_\pm = -\cos (\phi )$ (where again 
$\boldsymbol {v}_\pm$ is the unit outward normal to the upper and lower channel walls), then defining 
$\partial _{{\textit {wall}}}$ as the boundary of the wall, the wall energy is
\begin{equation} E_{{wall}} = \oint_{ \partial_{\textit{wall}}} \boldsymbol{w} \boldsymbol{\cdot} \,\mathrm{d} \boldsymbol{r}. \end{equation}
For the upper and lower walls described by 
$z = \pm 1/2 + B_{\pm }(y)$, we can choose, for example, 
$w_x = w_z = 0$ and 
$w_y = -x \cos \phi \sqrt {1 + B_\pm (y)^2}$. The integral around the closed curve is implemented internally by Surface Evolver along the edges specified by the user, with their orientation defined such that the unit normal is outward-pointing.
Boundary condition (2.4) is imposed by constraining the contact-line vertices to lie on the upper wall of the channel. This is a local condition on each vertex.
The fixed volume of liquid 
$V_L$ is handled in Surface Evolver as a global constraint on the possible energy configurations that the surface can take; that is, it removes one degree of freedom from the problem.
The mesh refinement is handled by Surface Evolver using a basic subdivision; we also equiangulate the mesh after each iteration. We converge to an energy minimum using the following process.
(i) Iterate on a fixed mesh until the solution is accurate to a specified tolerance.
(ii) Refine the mesh and check the difference between the energy on the new mesh and the old mesh. While the difference is greater than a specified tolerance, repeat step (i).
We use a tolerance of 
$10^{-6}$ for the accuracy of the solution on each mesh and the energy difference between meshes. We ensure that a global minimum has been reached by using a second-order gradient-descent method to check for positive eigenvalues near the equilibrium.
Appendix D. Numerical solution of the linear problem
We solve the linear problem (2.8)–(2.12) with Gaussian boundary data 
$b_\pm (y) = \pm \exp (-(y-y_c^\pm )^2/s^2)$ in a rectangular domain 
$-W \leq y \leq W$, 
$-{{\tilde {\theta }}} \leq \theta \leq {{\tilde {\theta }}}$. For general 
$y_c^\pm$, we integrate the Helmholtz equation (2.8) using second-order-accurate central finite differences with step lengths 
$\Delta y$ and 
$\Delta \theta$ in the 
$y$ and 
$\theta$ directions, respectively. We denote the value of the solution 
$f$ at 
$y = k \Delta y$, 
$\theta = j \Delta k$ by 
$f_k^j$ for 
$0 \leq k \leq M+1$, 
$0 \leq j \leq N+1$, so that 
$y = W$ is approximated by 
$(M+1)\,\Delta y$ and 
$\theta = {{\tilde {\theta }}}$ is approximated by 
$(N+1)\,\Delta \theta$. We discretise the Helmholtz equation (2.8) on the interior of the grid as
\begin{align} &\frac{1}{\Delta y^2}\,f_{k+1}^j + \frac{1}{\Delta y^2}\,f_{k-1}^j + \frac 1{\Delta \theta^2 R^2}\,f_k^{j+1} +\frac 1{\Delta \theta^2 R^2}\, f_k^{j-1} + \left( \frac{1}{R^2} - \frac{2}{\Delta \theta^2 R^2} - \frac{2}{\Delta y^2}\right)f_k^j = p, \nonumber\\ &\quad (1 \leq k \leq M, \ 1 \leq j \leq N). \end{align}We use the boundary conditions (2.9) and (2.10) to show that
$$\begin{gather} 2 \Delta \theta \sin(\theta_{N+1})\,f_k^{N+1} - 2 \Delta \theta\, b_+(y_k) +(3f^{N+1}_k - 4f^N_k +f^{N-1}_k) \cos(\theta_{N+1}) = 0, \end{gather}$$
$$\begin{gather}2 \Delta \theta \sin(\theta_{0})\,f^{0}_k + 2 \Delta \theta\,b_-(y_k) + (3f^{0}_k - 4f^1_k +f^{2}_k) \cos(\theta_{0}) = 0; \end{gather}$$meanwhile, the Neumann boundary conditions (2.11) give
\begin{equation} -3f_0^j + 4f_1^j -f_2^j = 0,\quad 3f_{M+1}^j - 4f_M^j +f_{M-1}^j = 0 \quad \mbox{for} \ 1 \leq j \leq N. \end{equation}
We then obtain a system of 
$(N+1)(M+1)$ equations that we solve subject to the volume constraint (2.12), which we discretise using Simpson's rule.
The discretised system is solved using a direct solver that takes advantage of the sparsity in the matrix structure. The grid size is chosen so that the solution to the discretised finite difference system at each grid point is accurate to three decimal places compared to the truncated analytical series solution, which is known at each grid point to high accuracy (truncated terms had size 
$O(10^{-16})$, see § 3.3).
Appendix E. The catenoid problem: governing equations and solution
Consider a catenoid with solid–liquid contact angle 
$0 \leq \phi < {\rm \pi}/2$ in a rectangular channel. Note that we use the term ‘catenoid’ here to describe the general shape of the interface as an ‘inverted droplet’; however, it need not be a surface of zero mean curvature. The catenoid interface is described by arc-angle coordinates 
$(r(t), z(t), \theta (t))$ for 
$-t_0 \leq t \leq t_0$ and is axisymmetric with respect to the azimuthal angle 
$\varphi$ in cylindrical polar coordinates 
$(r, \varphi, z)$. We take 
$t=0$ to be at 
$z=0$ as shown in figure 8 so that 
$(r(0), z(0), \theta (0) )= (r_0, 0, {\rm \pi}/2)$. Meanwhile, the channel walls are at 
$t = \pm t_0$ so that 
$(r(t_0), z(t_0), \theta (t_0)) = (R_d, 1/2, \phi )$ and 
$(r(-t_0), z(-t_0), \theta (-t_0)) = (R_d, -1/2, {\rm \pi}-\phi )$. The catenoid is symmetric about 
$z=0$, therefore without loss of generality we can consider the interface from 
$0\leq t \leq t_0$. Then 
$r$ and 
$z$ depend implicitly on 
$t$ as
$$\begin{gather} \frac{\mathrm{d} r}{\mathrm{d} t} = \cos \theta, \quad \frac{\mathrm{d} z}{\mathrm{d} t} = \sin\theta, \quad 0 \leq t \leq t_0, \end{gather}$$
$$\begin{gather}r(0) = r_0, \ r(t_0) = R_d, \quad z(0) = 0, \ z(t_0) = \tfrac{1}{2}. \end{gather}$$
The unit normal to the interface pointing into the vapour phase at 
$\varphi = {\rm const}.$ is given by 
$\hat {\boldsymbol {n}} = \sin \theta \, \hat {\boldsymbol {r}} -\cos \theta \,\hat {\boldsymbol {z}}$. Thus the Young–Laplace equation is
\begin{equation} \Delta p = \boldsymbol{\nabla} \boldsymbol{\cdot} {\hat{\boldsymbol{n}}} = \cos \theta\,\frac{\partial \theta}{\partial r} + \frac{\sin \theta}{ r} + \sin \theta\,\frac{\partial \theta}{\partial z} = \frac{\mathrm{d} \theta}{\mathrm{d} t} + \frac{\sin \theta}{ r}, \end{equation}so that the final equation in the system is
$$\begin{gather} \frac{\mathrm{d} \theta}{\mathrm{d} t} = \Delta p -\frac {\sin \theta}{r} , \quad 0 \leq t \leq t_0, \end{gather}$$
$$\begin{gather}\theta(0) = \frac {\rm \pi}2, \quad \theta(t_0) = \phi. \end{gather}$$
The ODEs (E1) and (E4), together with the boundary conditions, form a boundary-value problem for the arc-angle components 
$r, z, \theta$.
A catenoid in a rectangular channel with height 
$-\tfrac 1 2 \leq z \leq \tfrac 1 2$ and solid–liquid contact angle 
$\phi$.
A large-radius asymptotic solution (
$r_0\gg 1$) to this system can be found by writing
$$\begin{gather} r(t) = r_0 + r_1(t) + \frac{r_2(t)}{r_0} + \cdots, \quad \theta(t) = \theta_0(t) + \frac{\theta_1(t)}{r_0} + \cdots, \end{gather}$$
$$\begin{gather}z(t) = z_0(t) + \frac{z_1(t)}{r_0} + \cdots, \quad \Delta p = (\Delta p)_0 + \frac{(\Delta p)_1}{r_0} + \cdots, \end{gather}$$
$$\begin{gather}t_0 = L_0 + \frac 1 {r_0}\,L_1 + \cdots. \end{gather}$$
Solving the leading-order problem, we find from the leading-order approximation (
$\theta _0=(\Delta p)_0t+{\rm \pi} /2$, 
$z_0=(\sin (\Delta p)_0 t))/(\Delta p)_0$) that the pressure difference across the catenoid is
\begin{equation} (\Delta p)_0 ={-}2 \cos \phi, \end{equation}
which is consistent with the pressure difference of the unperturbed static liquid–vapour meniscus in the rectangular channel, while the leading-order approximation to the catenoid radius and the endpoint of the curve 
$t_0$ is
\begin{equation} r_1(t) = \frac{\cos((\Delta p)_0 t)}{(\Delta p)_0} - \frac 1 {(\Delta p)_0}, \quad L_0 = \frac{2 \phi - {\rm \pi}}{2}. \end{equation}
Solving the 
$O(r_0^{-1})$ problem, we find that
\begin{equation} (\Delta p)_1 ={\frac {\sin \tilde{\theta} \cos \tilde{\theta} +\tilde{\theta}}{2\sin \tilde{\theta} }}. \end{equation}
Thus, eliminating 
$r_0$ from truncated expansions for 
$r(t_0)$ and 
$\Delta p$, the expression
\begin{equation} R_d \approx \frac{(\Delta p)_1}{\Delta p - (\Delta p)_0 } + r_1(L_0) \end{equation}
gives the large-radius approximation for the catenoid for any given 
$\Delta p$.
We can also solve for catenoids with smaller radii numerically. First, we eliminate 
$t$ to obtain a nonlinear boundary-value problem where the unknown radius 
$R_d$ is to be determined as part of the solution for given 
$\Delta p$ and 
$\phi$:
$$\begin{gather} \frac{\mathrm{d} \theta}{\mathrm{d} z} = \frac{\Delta p}{\sin \theta} - \frac{1}{ r}, \quad \frac{\mathrm{d} r}{\mathrm{d} z} = \frac{\cos \theta}{\sin \theta}, \end{gather}$$
$$\begin{gather}r \left(\tfrac 1 2 \right) = R_d, \quad \theta(0 ) = \frac {\rm \pi}2, \quad \theta\left(\tfrac 1 2 \right) = \phi. \end{gather}$$
We solve this problem numerically using the Matlab routine ‘bvp4c’ (Kierzenka & Shampine Reference Kierzenka and Shampine2001), thus for any static meniscus with a given pressure difference (mean curvature) and contact angle, we can find the radius 
$R_d$ of the circular contact line of the catenoid that has the same pressure difference (mean curvature). The relationship between pressure difference and catenoid radius is shown in figure 9. It matches closely the asymptote (E12) for 
$R_d\gg 1$.
The maximum radius 
$R_d$ of a catenoid for varying pressure difference 
$\Delta p$, for a contact angle 
$\phi = 45^\circ$. The solid line denotes the numerical solution for 
$R_d$. The dashed line denotes the asymptotic solution (E12). The location of the asymptote is at 
$\Delta p = -2\cos ({\rm \pi} /4) \approx -1.414$.
The large-radius catenoid solution describes the curved meniscus shape (4.1) far from the wall perturbation, as we now demonstrate. The circular contact line is described locally by a parabola. To see this, take a catenoid with a contact line of radius 
$R_d$ centred on 
$x=x_0$, 
$y=W$. Its contact line lies along 
$(x-x_0)^2+(y-W)^2=R_d^2$. Because the solution is translationally invariant, 
$x_0$ may be chosen such that the catenoid passes through 
$x=0$, 
$y=y_c$. When 
$W\ll R_d$ and the centre of the catenoid lies in 
$x>0$, we may describe the base of the contact line using
\begin{gather} x=x_0-\sqrt{R_d^2-(W-y)^2} \end{gather}
\begin{gather} \approx x_0 - R_d\left(1-\frac{(W-y)^2}{2R_d^2}+\cdots \right) \approx C_0 + \frac{1}{2R_d}(W-y)^2, \end{gather}
with 
$C_0$ a constant. Therefore the contact-line displacement can be matched to (4.1) by choosing
\begin{equation} \epsilon p = \frac{\cos \tilde{\theta} \sin\tilde{\theta}+\tilde{\theta}}{2 R_d \sin\tilde{\theta}}=\frac{(\Delta p)_1}{R_d}, \end{equation}
and therefore approximates the contact line in 
$y_c< y\leq W$ when 
$\vert y-y_c \vert \gg s$. This pressure–radius relationship is consistent with the leading-order relationship found via the asymptotic expansion (E11), with 
$r_0 \approx R_d$.
We can then assess the limit 
$W\sim \epsilon ^{-1}$, for which the linearisation approximation of § 2.1 formally breaks down. Recall that the contact-line displacement 
$\epsilon x_p^\pm$ is 
$O(\epsilon pW^2)$ with 
$p=O(1/W)$ (from (3.1)), so that the contact-line displacement is 
$O(1)$, and 
$p=O(\epsilon )$ for 
$W\sim \epsilon ^{-1}$. However, the contact line retains a radius of curvature that is large compared to 
$W$ (
$R_d=O(1/\epsilon ^2)$, from (E12) with 
$\Delta p - (\Delta p)_0=\epsilon p$), allowing the parabolic approximation (E14) to be used. Thus the parabolic description (4.1) remains appropriate in this limit (figure 4), because of the structure of the catenoid solution. In contrast, larger-amplitude wall perturbations will cause 
$R_d$ to fall towards the size of 
$W$, pushing the contact line towards a more circular shape away from perturbations.
Appendix F. Sharp ridges and grooves
It is well known that a sharp wedge or groove can drive large contact-line displacements (Concus & Finn Reference Concus and Finn1969), and so far we have restricted attention to Gaussian perturbations (3.2) having width 
$s$ no smaller than 
$O(\epsilon ^{1/2})$, where 
$\epsilon$ is the wall-perturbation amplitude.
We now examine empirically what happens to the contact-line solution for the Gaussian perturbations 
$b_\pm (y) = \exp (-y^2/s^2)$ with 
$s \to 0$. Figure 10(a) shows the upper contact-line solutions for a channel-volume-preserving perturbation with decreasing 
$s^2$, with the narrowest computed perturbation having 
$s^2=0.0025$ in a channel of half-width 
$W=2$. Linear solutions computed using the series solution (3.3) are compared to Surface Evolver solutions, which are converged to an accuracy of 
$10^{-8}$ using the process outlined in Appendix C. The amplitude of the contact-line displacement increases as the perturbations become narrower; the linear model underpredicts the nonlinear Surface Evolver solution, indicating that nonlinear effects become important as the perturbations become sharper, particularly once 
$s^2$ approaches 
$\epsilon =0.01$. We also note that the far-field solution does not quite have zero curvature (as the linear model predicts for a channel-volume-preserving perturbation); this again is a nonlinear effect. Plotting the maximum displacement of the contact line at 
$y=0$ (figure 10b) shows that the amplitude of the displacement scales like 
$\log (1/s)$, suggesting that blowup may be possible even for analytic boundary forcing.
(a) The upper contact-line displacement 
$\hat {x}_p^+(y)$ for a channel-volume-preserving Gaussian perturbation 
$B_\pm (y)=0.01\exp (-y^2/s^2)$, with 
$s^2$ taking values 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1, 0.25 and 0.5. The black lines denote the linear solution, computed using the series solution (3.3), while the coloured dots denote the Surface Evolver solution. The contact angle is 
$\phi = 85^\circ$, and the channel half-width is 
$W=2$. (b) The upper contact-line displacement at 
$y=0$, 
$\hat {x}_p^+(0)$, with a logarithmic scale on the 
$x$ axis (values of 
$s$). The crosses denote the values of 
$\hat {x}_p^+(0)$ computed using the Surface Evolver solution; the circles are the corresponding values for the linear solution. The dashed line is 
$x_p^+(0)=0.0128 \log s-0.0057$.
A more extreme response can be expected for smaller contact angles and less smooth forcing. Consider the case in which the ridge or groove has small amplitude and narrower width (not necessarily Gaussian, but effectively satisfying 
$s^2\ll \epsilon$). More specifically, setting 
$s$ to zero, suppose that the lower wall shape (
$b_-$), say, has a discontinuity in a derivative at 
$y=0$, such that 
$b_-(y)=0$ for 
$y<0$, and 
$b_-(y)=y^\gamma$ for 
$y>0$. Then 
$\gamma =0$ corresponds to a step in 
$b_-$, 
$\gamma =1$ corresponds to a corner (a discontinuity in slope) and 
$\gamma =2$ corresponds to a jump in the curvature of the wall. The linearised curvature of the nearby gas–liquid interface, described in general by the Helmholtz equation (2.8), can be expected to be approximated in the neighbourhood of the discontinuity by 
$\nabla ^2 f\approx 0$. In the fully wetting case, for example, with 
${\tilde {\theta }}={\rm \pi} /2$, 
$b_-$ imposes 
$f$ along the wall (
$\theta =-{\rm \pi} /2$) via (2.9), and the wall normal derivative 
$f_n$ determines the contact-line displacement via (2.13a,b). Introduce polar coordinates 
$(\varrho,\vartheta )$ centred on 
$y=0$, such that 
$\vartheta =0$ (
$\vartheta ={\rm \pi}$) lies along the wall for 
$y>0$ (
$y<0$), and consider first the case of a step (
$\gamma =0$). Then 
$f(\varrho,\vartheta )=-(1-\vartheta /{\rm \pi} )$ provides a local solution to Laplace's equation subject to the forcing condition 
$b_-(y)=H(y)$, where 
$H$ is a Heaviside function. The corresponding wall normal derivative 
$f_n$ is then proportional to 
$1/y$, indicating that the contact line will be displaced in opposite directions on either side of a step. Further cases follow by integrating with respect to 
$y$, so that 
$f_n\propto \log \vert y\vert$ for 
$\gamma =1$ (the contact line will be displaced along the axis of a corner) and 
$f_n\propto y\log \vert y\vert -y$ for 
$\gamma =2$. These approximate solutions suggest that a very sharp step or a corner in wall shape, even if smoothed over a very short length scale 
$s$, will cause substantial deflection of the contact line (violating the linearisation approximation), while a jump in wall curvature will bend the contact line sufficiently for it to have infinite slope with respect to 
$y$, while remaining continuous.
In summary, and as indicated by figure 10, nonlinear effects will have a leading-order role close to the ridge or groove whenever the wall shape is sufficiently sharp.
Open access
