1. Introduction
The transport and manipulation of a thin liquid film arises in a wide variety of technological applications, such as coating and printing (Kalliadasis, Bielarz & Homsy Reference Kalliadasis, Bielarz and Homsy2000; Wierschem, Scholle & Aksel Reference Wierschem, Scholle and Aksel2002; Decré & Baret Reference Decré and Baret2003; Weinstein & Ruschak Reference Weinstein and Ruschak2004), drying (Cairncross, Francis & Scriven Reference Cairncross, Francis and Scriven1996) and heat exchange processes (Das, Choi & Patel Reference Das, Choi and Patel2006). In the biological context, examples of thin-film dynamics can be found in the membranes of mammalian lungs (Grotberg Reference Grotberg1994; Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997) and tear films in the eye (Braun Reference Braun2012).
Since the formulation of the analytical Nusselt solution (Nusselt Reference Nusselt1916), the transport of thin films driven by gravity has been extensively studied, with the focus on the nonlinear interaction between the fluid flow and the surface wave structure (Kapitza Reference Kapitza1948; Kapitza & Kapitza Reference Kapitza and Kapitza1949; Benjamin Reference Benjamin1957; Yih Reference Yih1963; Floryan, Davis & Kelly Reference Floryan, Davis and Kelly1987; Pozrikidis Reference Pozrikidis1988; Prokopiou, Cheng & Chang Reference Prokopiou, Cheng and Chang1991; Chang Reference Chang1994; Liu, Schneider & Gollub Reference Liu, Schneider and Gollub1995; Kalliadasis et al. Reference Kalliadasis, Bielarz and Homsy2000; Samanta Reference Samanta2022). By contrast, thin films driven by a peristaltic surface have begun to garner attention relatively recently, inspired by biological systems. For instance, Joo et al. (Reference Joo, Jung, Lee, Cowie and Takagi2020) uncovered an unusual feeding strategy of apple snails that utilise muscular undulations of their feet to drive thin, free surface flows and collect floating food particles. Motivated by this observation, Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023) built a simple robotic undulator that generates travelling waves to pump the liquid film. Their experiments have demonstrated that the pump rate of the thin film varies non-monotonically with the speed of the undulator's travelling wave. They also built a lubrication model that successfully rationalises their experimental observations.
Despite the success of their model, the work by Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023) solely focuses on the regime in which viscous effects dominate over inertia. However, as illustrated in gravity-driven thin films, finite inertia can drastically alter the dynamics of free surface flows, pointing to a new physical regime of interest. For instance, Malamataris, Vlachogiannis & Bontozoglou (Reference Malamataris, Vlachogiannis and Bontozoglou2002) questioned the validity of the parabolic velocity profile for a flow over an inclined plane as inertia becomes non-negligible. Their analysis demonstrated that, with finite inertia, large solitary waves could develop on the free surface and also lead to the occurrence of a counterflow in the region of minimum fluid thickness. This theoretical finding was later corroborated experimentally by Tihon et al. (Reference Tihon, Tovchigrechko, Sobolik and Wein2003).
In the present manuscript, we extend the study of Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023) to include the effects of inertia by conducting new experiments with a robotic undulator. In the experiments with non-negligible inertial effects, the pump rates show a clear deviation from the viscous-dominated regime, which cannot be explained by the original lubrication model. To explain the experimental observations, we develop a new theoretical model for the thin-film flow. To account for inertia in a mathematically tractable manner, we employ an asymptotic expansion on the velocity and pressure fields (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville1998, Reference Ruyer-Quil and Manneville2000; Scheid, Ruyer-Quil & Manneville Reference Scheid, Ruyer-Quil and Manneville2006), under the condition of a Reynolds number around 
$O(1)$. Similar asymptotic expansions have been previously utilised to study the effects of inertia and the effects of non-uniform surface on gravity-driven thin films (Benney Reference Benney1966; Gjevik Reference Gjevik1970; Nakaya Reference Nakaya1975; Wang Reference Wang1981). The method of expansion allows us to linearise the inertial terms in the Naiver–Stokes equation, so that we can find the analytical expressions for the velocity field and subsequently solve for fluid fluxes and free surface profiles numerically. Furthermore, we perform an additional asymptotic expansion on the free surface shape in the limit of small interfacial deformations, which yields analytical solutions of the fluid flux and gives us more precise insight into the role of inertia on free surface pumping.
Overall, our new theoretical model successfully reproduces the pump rates from the experiments with no fitting parameters, in the regimes with finite and negligible inertia. With our simplified model, we are able to more systematically probe the role of inertia, gravity and viscosity on the free surface shape, going beyond the experimental data. The analysis is able to elucidate the dependence of the free surface wave structure on the relevant non-dimensionalised numbers as well as establish the connection between the pump rate and the free surface shape. The paper is organised as follows: we introduce the experimental set-up in § 2, while the formulation of the mathematical model and the numerical method are discussed in § 3. Section 4 comprises the results of the numerical simulations in comparison with the experimental observations. In § 4.3, the analytical solutions of the fluid fluxes and the free surface profiles are provided to qualitatively rationalise the relation between free surface structures and the fluid fluxes. Finally, we conclude the paper with the summary and discussion in § 5.
2. Experiments
To create a surface with wave-like movements, we engineer an initial mechanical system depicted in figure 1(a). This system is identical to that of Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023) which employs a helix's rotational motion to induce a sinusoidal travelling wave pattern on a connected surface such that the wave speed 
$V_w=\omega \lambda$. The helix features a constant pitch and radius which fixes the wavelength (
$\lambda$) and amplitude (
$\delta$) of the travelling wave; the helix's rate of rotation is given by 
$\omega$. The helix is driven by a low-power servo motor which controls 
$\omega$. The helix is embedded within multiple rectangular slits which are connected by a thin (
$\sim$0.1 mm) membrane. Rotation of the helix gives rise to a rhythmic undulation of this membrane forming the travelling wave. All components of this device were 3D printed with 
$\lambda =50$ mm and 
$\delta =2.5$ mm. On a larger length scale, an alternate method of generating a pedalling-like motion of the solid surface was demonstrated by Vivanco et al. (Reference Vivanco, Cartwright, Ledesma Araujo, Gordillo and Marin2021) who divided the bottom boundary into small elements with prescribed individual motion in an elliptical orbit.
(a-i) The schematic of the robotic undulator consisting of a motor that rotates a helix encased inside a thin membrane. (a-ii) The PIV measurement shown with the free surface on top and the undulating solid boundary on the bottom. (a-iii) The undulator profile over one period of undulation. (b) Instantaneous flow rates 
$Q_i$ as a function of time 
$t$: for silicone oil (blue) at 
$V_w = 90.4\ {\rm mm}\ {\rm s}^{-1}$ and 
$H =6.8$ mm and for water–glycerine (orange) at 
$V_w = 78.54\ {\rm mm}\ {\rm s}^{-1}$ and 
$H = 6.3$ mm. (c) Time-averaged flow rate, 
$\langle Q \rangle$ is plotted against the wave speed 
$V_w$ scaled with 
$\epsilon ^2 H$ for silicone oil (blue) at 
$H = 6.8$ mm and water–glycerine (orange) at 
$H = 6.3$ mm. The slope of the dashed line is 1.5.
The undulator is affixed to the bottom of a fluid-filled tank containing silicone oils (density 
$\rho =970 \ {\rm kg}\ {\rm m}^{-3}$, viscosity 
$\mu =0.97 \ {\rm Pa}\ {\rm s}$, interfacial tension 
$\sigma =0.021 \ {\rm N} {\rm m}^{-1}$) and glycerine–water mixtures (
$\rho =1164 \ {\rm kg}\ {\rm m}^{-3}$, 
$\mu =0.133 \ {\rm Pa}\ {\rm s}$, 
$\sigma =0.067 \ {\rm N} {\rm m}^{-1}$). The depth of the liquids in the tank is maintained to such levels that only a thin layer (
$H$) remains on top of the undulator (cf. figure 1a). Across all our experiments 
$H$ varies between 4.3 and 8 mm which are much smaller than the wavelength of the undulations, 
$\lambda$. In this work, we will specifically focus on the experiments with silicone oil at 
$H = 6.8$ mm and water–glycerine at 
$H = 6.3$ mm. Travelling waves on the undulator create a directional transport of liquid within the thin film. To minimise boundary effects, we employ a spacious, acrylic tank measuring 
$61\ {\rm cm} \times 46\ {\rm cm}$.
We focus on characterising the flow rate within the thin film by performing particle image velocimetry (PIV) analysis. To this end, we introduce hollow glass particles of diameter 10 
${\rm \mu}$m and illuminate these tracer particles with a 520 nm 1 Watt laser sheet. The continuous laser sheet is positioned above the liquid layer to illuminate a longitudinal plane that passes through the centre of the undulator. Then, we visualise the motion of the illuminated particles with a high-speed camera. A typical velocity field is shown in figure 1(a-ii), where the colour coding is based on the horizontal component of the velocity; blue is the portion of liquid with the horizontal velocity in the same direction as 
$V_w$ whereas liquid in the red region has the velocity opposite to 
$V_w$.
To estimate the net transport of liquid by an undulator, we integrate the horizontal component of the fluid velocity at the mid-point of the undulator across the film thickness, which yields an instantaneous flow rate 
$Q_i(t)$. Figure 1(b) shows the plot of 
$Q_i(t)$ as a function of time in silicon oil (blue) at 
$V_w = 90.4\ {\rm mm}\ {\rm s}^{-1}$ and in a water–glycerine mixture (orange) at 
$V_w = 78.54\ {\rm mm}\ {\rm s}^{-1}$. The instantaneous flow rate fluctuates with the same time frequency as the undulator but has a positive offset. This positive offset confirms the ability of the undulator to maintain a net flow rate over a period of oscillation. Thus, we introduce a time-averaged flow rate as
\begin{equation} \langle Q \rangle =\frac{1}{\tau}\int_0^{\tau}Q_i(t)\,\mathrm{d}t, \end{equation}
where the time scale 
$\tau$ is chosen to be larger than one period to avoid any end effect of sinusoidal instantaneous flow rates. Figure 1(c) shows the flow rate vs 
$\epsilon ^2V_w H$ measured in two different fluids; silicone oil (blue) at 
$H = 6.8$ mm and a water–glycerine mixture (orange) at 
$H = 6.3$ mm. The role of the geometric factor 
$\epsilon \equiv \delta /H$ in this plot will be clarified in the subsequent sections.
At lower values of 
$\epsilon ^2V_w H$ (i.e. 
$\lesssim 50$), 
$\langle Q \rangle$ is shown to increase linearly with the wave speed for both the glycerine–water mixture and silicone oil. This liquid-independent regime is a consequence of the negligible deformation of the free surface at smaller 
$V_w$ and follows the relation 
$\langle Q \rangle =(3/2)\epsilon ^2V_w H$, as previously described by Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023). At higher values of 
$\epsilon ^2V_w H$, qualitatively different behaviours are observed for the two liquids; while the pump rate of the silicone oil gradually increases with the wave speed, the glycerine–water mixture exhibits a non-monotonic flow rate. In the rest of the paper, we present a theoretical model that captures these different behaviours.
3. Theory
We construct a mathematical model based on the experiments. In the two-dimensional system illustrated in figure 2(a), we consider the flow in a film of a Newtonian liquid on an undulating boundary that undergoes a periodic deformation. The profiles of the undulating boundary and the free surface are denoted as 
$h_1(X,t)$ and 
$h_2(X,t)$, respectively, so that 
$h_2-h_1$ corresponds to the thickness of the thin film. For simplicity, we model the shape of the undulating surface as 
$h_1(X,t)=\delta \sin [2{\rm \pi} (X-V_wt)/\lambda ]$, which describes the time-dependent periodic deformation with the amplitude 
$\delta$, the wavelength 
$\lambda$ and the speed of the travelling wave 
$V_w$. Consistent with the experimental set-up, the 
$X$-coordinate is defined positive in the direction of the travelling wave. We acknowledge that the actual shape generated by the robotic undulator does not follow a simple sine wave, as shown in figure 1(a) iii and also previously illustrated by Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023).
Schematic of a thin-film flow driven by an undulating surface with wave speed 
$V_w$ and wavelength 
$\lambda$ (a) in the laboratory frame and (b) in the wave frame. Here, the shape of solid boundary is denoted as 
$h_1$, while 
$h_2$ corresponds to the free surface profile.
To eliminate the time-dependent effects, we consider a frame moving with the wave speed 
$V_w$ (Chan, Balmforth & Hosoi Reference Chan, Balmforth and Hosoi2005; Lee et al. Reference Lee, Bush, Hosoi and Lauga2008). For the steady-state assumption to be valid, the deformation of the fluid–fluid interface and corresponding fluid flow must also move with 
$V_w$, which is qualitatively supported by the fact that the oscillation frequency of 
$Q_i$ in figure 1(b) matches that of the undulator. We will further comment on the validity and limitation of our steady-state assumption later in this section. We use a Galilean transformation that relates the laboratory coordinates 
$(X, Y)$ to the wave-frame coordinates 
$(x, y)$, such that
\begin{equation} x =X - V_wt, \quad y = Y. \end{equation}Thus, the transformation of the velocity field can be written as
\begin{equation} u(x,y)=U(X,Y,t) - V_w, \quad v(x,y)=V(X,Y,t). \end{equation}
Here, 
$(u, v)$ denotes the wave-frame velocity field in the 
$(x, y)$ directions, whereas 
$(U,V)$ is the laboratory-frame velocity field in 
$(X,Y)$ (Shapiro, Jaffrin & Weinberg Reference Shapiro, Jaffrin and Weinberg1969). Note that we first solve for the flow system in the wave frame and then convert the results back to the laboratory frame to compare with the experimental observations.
Next, we non-dimensionalise the governing equations based on the following set of characteristic scales:
\begin{gather} (x,y)=\lambda(x^*,ay^*), \quad (h_1,h_2)=H(h_1^*,h_2^*), \end{gather}
\begin{gather}(u,v)={V_w}(u^*,av^*), \quad p= p^*\mu V_w\lambda/H^2, \end{gather}
where the asterisk denotes dimensionless variables, 
$p$ is the gauge pressure and 
$\mu$ is the liquid viscosity. Note that 
$a\equiv H/\lambda$, where 
$H$ is the mean thickness of the liquid film. In addition, the dimensionless shape of the undulating surface corresponds to 
$h_1^*(x^*) = \epsilon \sin {(2{\rm \pi} x^*)}$, where 
$\epsilon \equiv \delta /H$. Then, the linear momentum and mass conservation equations for the fluid flow can be expressed as
\begin{gather} \frac{\partial u^*}{\partial x^*}+\frac{\partial v^*}{\partial y^*}=0, \end{gather}
\begin{gather}a{\widetilde{{Re}}}\left(u^* \frac{\partial u^*}{\partial x^*}+v^* \frac{\partial u^*}{\partial y^*}\right)={-}\frac{\partial p^*}{\partial x^*}+a^2\frac{\partial^2 u^*}{\partial {x^*}^2}+\frac{\partial^2 u^*}{\partial {y^*}^2}, \end{gather}
\begin{gather}a^3{\widetilde{{Re}}}\left(u^* \frac{\partial v^*}{\partial x^*}+v^* \frac{\partial v^*}{\partial y^*}\right)={-}\frac{\partial p^*}{\partial y^*}+a^2\left(a^2\frac{\partial^2 v^*}{\partial {x^*}^2}+\frac{\partial^2 v^*}{\partial {y^*}^2}\right)-\frac{Bo}{{{Ca} }}, \end{gather}
with Reynolds number 
${\widetilde {{Re}}}=a\rho V_w \lambda /\mu$, Bond number 
${{Bo}}=\rho g \lambda ^2/\sigma$ and the capillary number 
${{{Ca}}}=\mu V_w/a^3\sigma$. Here, 
$g$ is the gravitational acceleration. The 
$Bo/Ca$ term in (3.7) corresponds to the dimensionless body force on the fluid.
The boundary conditions in the wave frame can be expressed as
\begin{gather} u^*+2{\rm \pi}\epsilon a^2\cos(2{\rm \pi} x^*)v^*={-}1 \quad \text{at } y^*=h_1^*, \end{gather}
\begin{gather}-2{\rm \pi}\epsilon\cos(2{\rm \pi} x^*)u^*+v^*=0\quad \text{at } y^*=h_1^*, \end{gather}
\begin{gather}\left[1+a^2\left(\frac{\partial h_2^*}{\partial x^*}\right)^2\right]\frac{\partial u^*}{\partial y^*}+2a^2\frac{\partial h_2^*}{\partial x^*}\left(\frac{\partial v^*}{\partial y^*}-\frac{\partial u^*}{\partial x^*}\right)+O(a^4) =0 \quad \text{at } y^*=h_2^*, \end{gather}
\begin{gather}p^*+a^2\left(\frac{\partial h_2^*}{\partial x^*}\frac{\partial u^*}{\partial y^*}-2\frac{\partial v^*}{\partial y^*}\right)+O(a^4)={-}\frac{{Ca}^{{-}1}\dfrac{\partial^2 h_2^*}{\partial {x^*}^2}}{ \left[1+a^2\left(\dfrac{\partial h_2^*}{\partial x^*}\right)^2\right]^{3/2}} \quad \text{at } y^*=h_2^*. \end{gather}
Here, (3.8) and (3.9) correspond to the no-slip and no-flux boundary conditions on the undulator surface, respectively. On the free surface, we impose the tangential stress balance via (3.10) and the jump in normal stress due to surface tension via (3.11). For simplicity, we drop the asterisk and use 
$\partial _x$ and 
$\partial _y$ to denote partial differentiation with respect to the 
$x$ and 
$y$ coordinates in the following discussion.
Motivated by the previous studies of thin-film flows that include inertial effects (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville1998, Reference Ruyer-Quil and Manneville2000; Scheid et al. Reference Scheid, Ruyer-Quil and Manneville2006), we take the limit of 
$a\ll 1$ and apply an asymptotic expansion of the unknowns (i.e. 
$u, v,$ and 
$p$) in terms of 
$a$ (Roy, Roberts & Simpson Reference Roy, Roberts and Simpson2002)
\begin{gather} u(x,y)=\sum_{j=0}^{\infty} a^j u_j(x,y), \end{gather}
\begin{gather}v(x,y)=\sum_{j=0}^{\infty} a^j v_j(x,y), \end{gather}
\begin{gather}p(x,y)=\sum_{j=0}^{\infty} a^j p_j(x,y). \end{gather}
Following the formulation of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville1998), we do not asymptotically expand 
$h_2$ in order of 
$a$, so that the free surface shape may include the coupled effects from all orders. In addition to 
$a\ll 1$, we assume that 
$\widetilde {{Re}}\sim O(1)$, which ensures that the leading inertial terms in the 
$x$-momentum equation show up at 
$O(a)$. Hence, our model approach is expected to fail at 
$\widetilde {{Re}}\sim O(a^{-1})$ or larger, as the linearisation of the inertial terms becomes problematic in that limit. Furthermore, at 
$\widetilde {{Re}}\gg 1$, we have experimentally observed the formation of free surface waves whose speed is distinct from the wave speed of the undulating surface. The disconnect between the two wave speeds means that our steady-state assumption in the wave frame is no longer valid. Therefore, we emphasise that the current model is valid for 
${\widetilde {{Re}}}\sim O(1)$ and 
$a\ll 1$.
The boundary conditions also need to be decomposed to match the governing equations at each order. Note that we assume 
$\epsilon$ to be comparable to 
$a$. The curvature term in the normal stress boundary condition can be expanded as 
$\partial _{xx}h_2 (1-3a^2(\partial _{x}h_2)^2/2+O(a^4))$ by applying the binomial expansion. However, since our main purpose for the asymptotic expansion is to include the inertial effects into the model, we neglect the higher-order terms in the normal stress boundary condition and simplify (3.11) to 
$p=-\partial _{xx}h_2/{{{Ca}}}$. We acknowledge that the 
$O(a^2)$ term may become important in the regime where the free surface deformations are significant, which mostly falls outside our current regime of interest. Hence, we reasonably neglect the higher-order terms of the normal stress boundary condition in the rest of the discussion. The boundary conditions can be decomposed into
\begin{gather} u_0={-}1, u_1=u_2=\cdots=0 \quad\text{at } y=h_1, \end{gather}
\begin{gather}v_0=v_1=v_2=\cdots=0 \quad\text{at } y=h_1, \end{gather}
\begin{gather}\partial_y u_0=\partial_y u_1=\partial_y u_2=\cdots=0 \quad\text{at } y=h_2, \end{gather}
\begin{gather}p_0={-}\partial_{xx}h_2/{{{Ca}}}, p_1=p_2=\cdots=0 \quad\text{at } y=h_2. \end{gather}
Then, at the zeroth order, we find the following governing equations for 
$u_0, v_0$ and 
$p_0$:
\begin{gather} \partial_x u_0+\partial_y v_0=0, \end{gather}
\begin{gather}0={-}\partial_x p_0+\partial_{yy}u_0, \end{gather}
\begin{gather}0={-}\partial_y p_0-{Bo}/{{{Ca}}}. \end{gather}
We note that 
$u_0(x,y)$ is a second degree polynomial, and the solution does not include the inertial terms.
Thus, to study the inertial effects on the free surface flow, we expand the unknowns in the next two orders. Then, at 
$O(a^1)$, we obtain
\begin{gather} \partial_x u_1+\partial_y v_1=0, \end{gather}
\begin{gather}{\widetilde{{Re}}}(u_0\partial_x u_0+v_0\partial_y u_0)={-}\partial_x p_1+\partial_{yy}u_1, \end{gather}
\begin{gather}0={-}\partial_y p_1, \end{gather}
which yield the solution of 
$u_1$ as a sixth degree polynomial of 
$y$. The equations at 
$O(a^2)$ correspond to
\begin{gather} \partial_xu_2+\partial_y v_2=0, \end{gather}
\begin{gather}{\widetilde{{Re}}}(u_0\partial_x u_1+u_1\partial_x u_0+v_0\partial_y u_1+v_1\partial_y u_0)={-}\partial_x p_2+\partial_{xx}u_0+\partial_{yy}u_2, \end{gather}
\begin{gather}0={-}\partial_y p_2+\partial_{yy}v_0, \end{gather}
resulting in 
$u_2$ as a tenth degree polynomial of 
$y$. Note that both 
$u_1$ and 
$u_2$ are explicit functions of 
$\widetilde {{Re}}$. Finally, to obtain the resulting liquid flux 
$q$ through the liquid layer, we integrate the velocity field in 
$y$, so that
\begin{equation} q=\int_{h_1}^{h_2}\sum_{j=0}^{n} u_j(x,y)\, {\rm d} y.\end{equation}
Here, we replace the infinite sum in the asymptotic expansion with the partial sum up to the order 
$n$. We select 
$n = 2$ to provide an appropriate physical description of inertial effects with a reasonable computational cost (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville1998, Reference Ruyer-Quil and Manneville2000). In the current wave-frame coordinates that are independent of time, conservation of mass ensures that 
$q$ is constant.
Notably, 
$q$ and 
$h_2$ are unknown a priori. Hence, we simultaneously solve for 
$q$ and 
$h_2$ by imposing additional boundary conditions. As previously noted, 
$q$ and 
$h_2$ are not expanded in orders of 
$a$. This enables us to numerically compute the free surface shape and fluid flux that include the effects from all orders. Based on the assumption that the undulating surface undergoes periodic deformations, we apply periodic boundary conditions for 
$h_2$ at 
$x=0$ and 
$x=1$. The periodic boundary conditions over each wavelength are: 
$h_2(0)=h_2(1)$, 
$\partial _x h_{2}(0)=\partial _x h_{2}(1)$, 
$\partial _{xx} h_{2}(0)=\partial _{xx}h_{2}(1)$, 
$\ldots$, 
$\partial _{x^{n+2}} h_{2}(0)=\partial _{x^{n+2}}h_{2}(1)$ where 
$n$ is the number of terms in the partial sum. In addition, we assume that the total fluid area 
$A$ is conserved (i.e. 
$A=\int _{0}^{1} (h_2-h_1) \,{{\rm d}\kern0.06em x}=1$).
Then, the method of solving 
$q$ and 
$h_2$ can be summarised as the following steps, with the details of the numerical scheme included in Appendix A:
(i) Make an initial guess for the flux
$q$.(ii) Find
$h_2(x)$ by solving (3.28) numerically with the periodic boundary conditions.(iii) Check the total area
$A_i=\int _{0}^{1} (h_2-h_1) \,{{\rm d}\kern0.06em x}$.(iv) Change the initial guess
$q$ until the root of the function $f(q)=A-A_i$ reaches zero.
The numerical solutions of 
$h_2$ and 
$q$ are presented in §§ 4.1 and 4.2. Furthermore, we obtain the analytical solution of 
$h_2$, by asymptotically expanding it in the limit of a ‘nearly flat’ surface. This allows us to derive an explicit relationship between the flux, the free surface shape and the key dimensionless parameters. The analytical results are presented in § 4.3.
4. Results
4.1. Comparison with experiments
4.1.1. Silicone oil
We first extract the free surface shapes of silicone oil (SO) at 
$V_w=0.0177 \ {\rm m}\ {\rm s}^{-1}$ and compare them with our theoretical prediction. Figure 3 shows the experimental oil–air interface in solid lines at two different times in the presence of an undulating boundary. Note that the experimental profiles are accompanied by the instantaneous PIV measurements in the laboratory frame. Figure 3 also includes the free surface shape and the undulating surface profile from our model in dashed lines. In order to compare our steady state solution with the experimental data at two different times, we have shifted our solution by an empirical offset of 
$\delta _{off}\approx 0.873$, such that 
$h_1=\sin (2{\rm \pi} x+0.873)$ to fit the experimental image at 
$t=0.48\ {\rm s}$. We do not expect the match between theory and experiments to be perfect, as the experimental solid boundary cannot be accurately described by a simple sine wave. We quantify the differences in 
$h_1$ and 
$h_2$ between theory and experiments by computing their 
$L^2$ norm, or 
$\|h_1 - h_{1\exp }\|$ and 
$\|h_2 - h_{2\exp }\|$, where the subscript ‘exp’ denotes the profiles from the experimental measurements. We find that 
$\|h_1 - h_{1\exp }\|$ and 
$\|h_2 - h_{2\exp }\|$ are of the same order of magnitude and are approximately 
$0.1$ over the entire domain. Despite this discrepancy, the results in figure 3 indicate the crest of the theoretical free surface profile is slightly ahead of that of the solid boundary, which qualitatively coincides with the experimental observations.
Results of PIV for silicone oil at 
$V_w= 0.0177 \ {\rm m}\ {\rm s}^{-1}$ at two different instances of time: (a) at 
$t=0$ and (b) 
$t=0.48 \ {\rm s}$. The colour coding inside the thin film represents the horizontal component of the laboratory-frame velocity field that is taken to be positive in the direction of 
$V_w$. The black and the blue dashed line are the theoretical predictions of the dimensional free surface shape and the dimensional profile of the undulating boundary with (a) 
$h_1=\sin (2{\rm \pi} x)$ for 
$x\approx 0\unicode{x2013}0.81$ and (b) 
$h_1=\sin (2{\rm \pi} x+0.873)$ for 
$x\approx 0\unicode{x2013}0.81$.
To further delve into the free surface shapes, figure 4(a) shows the theoretical prediction for the free surface shapes at 
$V_w= 0.02$, 
$0.04$, 
$0.06$ and 
$0.08 \ {\rm m}\ {\rm s}^{-1}$. As 
$V_w$ increases, the shape of the free surface gradually starts to resemble the undulating surface. The physical mechanism behind free surface deformations can be understood by comparing the dimensionless numbers in the flow system. Based on the governing equations and boundary conditions (i.e. (3.5) to (3.11)), we note that the flow system is controlled by 
$Re, Bo$ and 
$Ca$, where we have defined a modified Reynolds number 
${Re} \equiv a{\widetilde {Re}}$ to account for the prefactor 
$a$ in front of 
${\widetilde {Re}}$ in (3.6). Since the magnitude of 
${Re}$ is around 
$O(10^{-2})$ as 
$V_w\sim O(10^{-2}\ {\rm m}\ {\rm s}^{-1})$, the inertial effect is insignificant, and the flow system is governed by 
${Ca}$ and 
$Bo$ only. Specifically, 
${Ca}/{Bo}\equiv \mu V_w /(a^3 \rho g \lambda ^2)$ represents the competition between viscous and hydrostatic effects, which increases linearly with increasing 
$V_w$. Hence, for smaller 
$V_w$ (or smaller 
${Ca}/{Bo}$), a relatively flat liquid–air interface can be observed since the hydrostatic effects are dominant. By contrast, with increasing 
$V_w$ (i.e. larger 
${Ca}/{Bo}$), the free surface deformations tend to increase, which also coincides with an increase in 
$Ca$, or reduced surface tension effects.
(a) The theoretical prediction of the free surface shapes for varying 
$V_w$. The different solid lines (coloured from light to dark green) correspond to the different values of 
$V_w=0.02$, 
$0.04$, 
$0.06$ and 
$0.08 \ {\rm m} {\rm s}^{-1}$. (b) The plot of laboratory-frame flux from the experiments (black dot) and theoretical predictions (blue lines) of SO. The error bars account for the uncertainty associated with the experimental measurement of 
$H$. The solid line and dashed line correspond to theoretical predictions including inertial effects (
$Q$) and under lubrication approximations (
$Q_{lub}$), respectively.
Next, we discuss how the liquid flux varies with 
$V_w$. We propose two different theoretical predictions of the liquid flux and compare them with the experimental results. The first prediction includes inertial effects, which is equivalent to 
$q$ in (3.28). The second prediction does not include inertial effect and is denoted as 
$q_{lub}$. In order to obtain 
$q_{lub}$, we follow the same procedure in § 3 but only include the terms at 
$O(a^0)$. Since all the inertial terms are neglected, the solution is identical to the results under lubrication approximations. However, 
$q$ and 
$q_{lub}$ are wave-frame fluxes, whereas the experimental results are in the laboratory frame.
To compare the theoretical predictions with the experimental results, we transform the theoretical results of the wave-frame flux into the laboratory frame by using the relation 
$U(X,Y,t) = u(x,y) + 1$. We find 
$Q_i(X_0,t) = q + (h_2(X_0,t)-h_1(X_0,t))$ by integrating the velocity from 
$Y=h_1(X_0,t)$ to 
$Y=h_2(X_0,t)$. Here, 
$Q_i(X_0,t)$ is the instantaneous laboratory-frame flux at a horizontal location 
$X_0$, examples of which were previously plotted in figure 1(b) from the PIV data. We time average 
$Q_i$ based on
\begin{equation} \frac{1}{\tau}\int_0^{\tau} Q_i(X_0,t)\, {\rm d}t = Q(X_0)= q+(\overline{h_2}(X_0)-\overline{h_1}(X_0)).\end{equation}
Due to the periodic nature of flow system, we note that 
$\overline {h_2}(X_0)-\overline {h_1}(X_0)=1$. Since the right-hand side of (4.1) is a constant, the laboratory-frame flux 
$Q$ is independent of the horizontal location 
$X$. Thus, we arrive at the following equation:
\begin{equation} Q= q+1. \end{equation}
The above transformation is also applicable under the lubrication assumption, such that 
$Q_{lub}=q_{lub}+1$.
Figure 4(b) shows the experimental results of 
$Q$, plotted with the theoretical predictions of 
$Q$ (solid line) and 
$Q_{lub}$ (dashed line) for the range of 
$V_w=0.01\unicode{x2013}0.1 \ {\rm m}\ {\rm s}^{-1}$. As the inertial effects are negligible (i.e. 
${Re}\ll 1$), there are only small differences between 
$Q_{lub}$ and 
$Q$, and they both effectively capture the overall decrease in the liquid flux with an increase in 
${Ca}/{Bo}$. As 
$V_w$ increases, the viscous effects gradually dominate the gravitational effects, which implies that the fluid becomes easier to undulate and move with the bottom profile, yielding a decrease in 
$Q$. However, despite the qualitative match between theory and experiments, theoretical 
$Q$ tends to consistently over-predict the experimental values due to the simplified nature of our model. Specifically, the two-dimensional model cannot capture any leakage flow from the sides of the undulator, which may explain why the experimental measurements tend to be lower than the idealised model results. Finally, while the differences 
$Q_{lub}$ and 
$Q$ are small, 
$Q_{lub}$ is systematically lower than 
$Q$ as 
$V_w$ is increased. This implies that the inertial effects tend to enhance the flux, which will be further explored in the next section.
4.1.2. Glycerine–water mixture
Next, we consider the case of choosing the glycerine–water (GW) mixture as the working liquid at 
$V_w=0.023 \ {\rm m}\ {\rm s}^{-1}$. As shown in figure 5, the amplitude of the surface wave is relatively small and the free surface shape is flat for both experiments (solid lines) and theory (dashed lines). We set 
$h_1=\sin (2{\rm \pi} x+ 0.698)$ and 
$h_1=\sin (2{\rm \pi} x+ 2.09)$ to best match the experimental undulating boundaries in figures 5(a) and 5(b), respectively. Consistent with the results in figure 3, 
$\|h_1 - h_{1\exp }\|$ and 
$\|h_2 - h_{2\exp }\|$ from figure 5 are again comparable to each other and of 
$O(10^{-1})$. To understand the resultant free surface shape, we note that the magnitude of 
$Ca/Bo$ for the GW mixture is an order of magnitude smaller than the SO case in figure 3, while the magnitude of 
$Ca$ is two orders of magnitude smaller than the SO data. Hence, as the hydrostatic and surface tension effects tend to dominate in the present case, the fluid–fluid interface remains relatively flat.
Results of PIV for the glycerine–water mixture at 
$V_w= 0.023 \ {\rm m}\ {\rm s}^{-1}$: (a) at 
$t=0$ and (b) 
$t=0.96 {\rm s}$. The black and the blue dashed lines are the theoretical dimensional free surface shapes and the dimensional profiles of the undulating boundary, which are given by (a) 
$h_1=\sin (2{\rm \pi} x+ 0.698)$ for 
$x\approx 0\unicode{x2013}0.81$, and (b) 
$h_1=\sin (2{\rm \pi} x+ 2.09)$ for 
$x\approx 0\unicode{x2013}0.81$.
Furthermore, in the experiments with the GW mixture, the magnitude of 
$Re$ ranges from 
$O(10^{-1})$ to 
$O(1)$, which means that the inertial effects are no longer negligible. However, due to the dominance of 
$Bo$ and the flat free surface in figure 5, the role of inertial effects remains unclear. For the following discussion, we will systematically change the wave speed to see how inertial effects may affect the free surface profile. In figure 6(a), we present the theoretical predictions of free surface profiles with 
$V_w= 0.03$, 
$0.06$, 
$0.09$ and 
$0.12 \ {\rm m}\ {\rm s}^{-1}$, which show that the interface deforms more as 
$V_w$ increases. At the same time, the free surface gradually becomes more out of phase with the undulating surface, as more clearly illustrated in figure 6(c). The zoomed-in plot of 
$h_2$ in figure 6(c) shows that the minimum point of the free surface (denoted as a triangle) moves to the left towards the peak of the undulating substrate. This means that, distinct from the SO case, the free surface profile no longer tends to conform to the undulating surface, as the wave speed is increased.
(a) The theoretical prediction of free surface shapes, where the different solid lines (coloured from light to dark green) correspond to the different values of 
$V_w=0.03$, 
$0.06$, 
$0.09$ and 
$0.12 \ {\rm m}\ {\rm s}^{-1}$. (b) Free surface shapes under under lubrication approximation for 
$V_w=0.03$, 
$0.06$, 
$0.09$ and 
$0.12 \ {\rm m}\ {\rm s}^{-1}$ (coloured from light to dark green). (c) The zoom-in plot of the free surface profiles from the lubrication model in (a); the red triangle corresponds to the minimum point in 
$h_2$. (d) The zoom-in plot of the free surface profiles from the lubrication model in (b); the red triangle corresponds to the minimum point in 
$h_2$. (e) The plot of laboratory-frame flux from the experiments (black dot) and theoretical predictions (blue lines) of GW. The error bars account for the uncertainty associated with the experimental measurement of 
$H$. The solid line and dashed line correspond to theoretical predictions including inertial effects (
$Q$) and under lubrication approximations (
$Q_{lub}$), respectively.
To isolate the effects of inertia on this phenomenon, we compute and plot the free surface profiles under lubrication approximations in figure 6(b). Although the deformation of the free surface profile is insignificant compared with the SO case, it gradually changes to become more in phase with 
$h_1$. This is evident in the zoomed-in plot of figure 6(d) that shows that the minimum point in 
$h_2$ gradually moves to the right with increasing 
$V_w$, in direct contrast to the model results in figure 6(c). The above observation implies two things: (i) the inertial effects enhance the amplitude of free surface deformations, and (ii) inertia causes the free surface to fall out of phase with 
$h_1$.
Next, we discuss how the free surface profile influences the fluid flux. Figure 6(e) shows the experimental results of 
$Q$ and theoretical predictions for 
$Q$ and 
$Q_{lub}$ for the range of 
$V_w=0.01\unicode{x2013}0.1 \ {\rm m}\ {\rm s}^{-1}$. We clearly observe that only the solution with the inertial effects is able to correctly capture the increase in 
$Q$ with the 
$V_w$, while 
$Q_{lub}$ gradually decreases with 
$V_w$. Hence, these results once again validate the role of inertia to increase the fluid flux. To connect 
$Q$ to free surface shapes, we recall that inertia causes the free surface deformations to be out of phase with the lower boundary, which suggests the potential interaction between the free surface and undulating substrate. At higher Reynolds numbers, the free surface tends to form a dip towards the peak of the undulating surface. The fluid is then ‘squeezed’ by the bottom boundary and free surface that move towards each other, which, thereby, increases the liquid flux.
4.2. The effects of ${Re}$ and ${Ca}/{Bo}$
While our results so far suggest that inertia enhances the fluid flux, we cannot clearly decouple the effects of inertia from viscous and hydrostatic ones, as both 
${Re}$ and 
${Ca}/{Bo}$ increase with 
$V_w$. To address this problem, we construct a phase map of the flux 
$Q$ by varying 
${Re}$ and 
${Ca}/{Bo}$ independently. We first set 
$a=0.126$, 
$H=6.3 \ {\rm mm}$ and 
$\sigma =0.02 \ {\rm N} \ {\rm m}^{-1}$ for all the simulations herein. While 
$\sigma$ drops out of 
${Ca}/{Bo}$, the value of 
$\sigma$ is needed to satisfy the normal stress boundary condition at the free surface (i.e. (3.11)). In the regime with large Bo (i.e. 
${Bo}\sim O(10^2)\unicode{x2013}O(10^3)$), we find that the effects of varying 
$\sigma$ on 
$Q$ are minimal and can be neglected.
To obtain 
$Q$ at given 
$Ca/Bo$ and 
$Re$, we start the simulation at 
$(Ca/Bo, Re)=(0.01,0)$ and gradually increase 
${Ca}/{Bo}$ (with an increment of 0.1) and 
${Re}$ (with an increment of 0.1). Once the values of 
$Q$ are obtained for the given range of 
${Ca}/{Bo}$ and 
${Re}$, we apply a Matlab function, griddata with the cubic method to interpolate our discrete simulation data and generate a smooth phase map. Notably, because 
$Re$ and 
${Ca}/{Bo}$ are both linearly dependent on 
$V_w$, their ratio is independent of 
$V_w$ and is given by
\begin{equation} \frac{{Re}}{{Ca}/{Bo}} = a^2H^3 g\left(\frac{\rho}{\mu}\right)^2, \end{equation}
where 
$a$, 
$H$ and 
$g$ are not fluid-dependent and remain fixed in all simulations. This implies that, once a working fluid is selected, 
$Q$ for the given fluid can be represented by a straight line passing through 
$(0,0)$ with a fixed slope, where 
${Re}/({Ca}/{Bo})$ corresponds to the slope.
The resulting 
${Ca}/{Bo}$–
${Re}$ phase maps of 
$Q$ are shown in figure 7; the phase map in figure 7(b) reproduces a smaller range of 
$Re$ from figure 7(a) for clarity. In addition, the two arrows in figure 7 correspond to the simulation results for the SO and the GW mixture, respectively. The phase maps in figure 7 reveal that 
$Q$ decreases with 
${Ca}/{Bo}$ for small 
$Re$, qualitatively similar to the SO case. However, for large 
$Re$, this trend is reversed, as 
$Q$ is shown to increase with 
${Ca}/{Bo}$, with the maximum 
$Q$ corresponding to the top-right corner of figure 7(a).
(a) Phase diagrams summarising the magnitude of 
$Q$ with varying 
${Re}$ and 
${Ca}/{Bo}$. The two arrows are identified as the fluid types SO and GW. The direction of the arrows corresponds to an increase in 
$V_w$. (b) A the zoomed-in plot of (a) that highlights the range of 
$Re=0\unicode{x2013}0.5$. The grey contour lines correspond to constant 
$Q$.
For an intermediate range for 
$Re$ (i.e. 
${Re}=0.1\unicode{x2013}0.5$), we can clearly see in figure 7(b) that 
$Q$ varies non-monotonically with 
$Ca/Bo$. To investigate this further, we plot 
$Q$ as a function of 
${Ca}/{Bo}$ for different values of 
$Re$ in figure 8. Figure 8(a) demonstrates that 
$Q$ monotonically decreases with 
${Ca}/{Bo}$ at small 
${Re}$ (see blue lines) and gradually exhibits a non-monotonic behaviour with 
${Ca}/{Bo}$ as 
${Re}$ is increased (see green lines). Then for 
${Re} > 0.5$, 
$Q$ is shown to increase monotonically with 
${Ca}/{Bo}$ in figure 8(b). Hence, the effects of increasing 
${Ca}/{Bo}$ on 
$Q$ are strongly dependent on 
$Re$, pointing to the potential mixed effects of 
${Ca}/{Bo}$ and 
$Re$ on 
$Q$.
(a) The theoretical predictions of 
$Q$ plotted against 
${Ca}/{Bo}$ with different values of 
${Re}$ from 0 to 0.5, and (b) from 0.6 to 0.8. (c) The theoretical predictions of 
$Q$ plotted against 
$Re$ with different values of 
$Ca/Bo$ from 0.2 to 0.8. The inset shows the predictions of 
$Q$ as a function of 
$Re$ for 
${Ca/Bo}=0.01$ (blue) and 
$0.1$ (green). The 
$y$-axis of the inset ranges from 
$Q=0$ to 0.3, while the range of 
$Re$ on the 
$x$-axis of the inset is the same as the main plot.
On the other hand, the phase diagrams in figure 7 appear to validate our earlier observation that inertia enhances the flux. To further corroborate this, we plot 
$Q$ as a function of 
$Re$ at constant 
$Ca/Bo$ in figure 8(c). The main panel in figure 8(c) shows a monotonic increase in 
$Q$ with increasing 
$Re$ for 
$Ca/Bo$ ranging from 0.2 to 0.8. This demonstrates positive effects of inertia on the flux. Surprisingly, at low 
$Ca/Bo$ (i.e. 
${Ca}/{Bo} = 0.01, 0.1$), the inset in figure 8(c) reveals that 
$Q$ tends to decrease as 
$Re$ approaches 
$O(1)$; this is also evident in figure 8(b) that shows smaller 
$Q$ for larger 
$Re$ near 
${Ca/Bo}\approx 0$. This unexpected reduction in 
$Q$ at relatively large 
$Re$ at low 
$Ca/Bo$ has two contributions. First, low 
$Ca/Bo$ (i.e. large gravitational effects) ensures minimal deformation of the fluid–fluid interface, which suppresses changes in 
$Q$ from the free surface shape. Second, inertia-induced changes in the velocity profile can enhance viscous dissipation from the higher-order terms, which tend to reduce 
$Q$ at low 
$Ca/Bo$. Apart from this small 
$Ca/Bo$ and large 
$Re$ regime, increasing 
$Re$ clearly enhances the overall flux. Hence, we will primarily focus on the non-monotonic behaviour of 
$Q$ with 
$Ca/Bo$ for non-zero 
$Re$ in the next section.
Finally, the results with the SO and GW mixture in § 4.1 have demonstrated the connection between the free surface shape and the resultant flux 
$Q$. Namely, in the limit of 
${Re}\ll 1$, increasing 
$Ca/Bo$ leads to the transition of the fluid–fluid interface from flat to that conforming to the solid surface, which yields a decrease in 
$Q$. On the other hand, for 
$Re \sim O(0.1)\unicode{x2013}O(1)$, increasing 
$Re$ (and 
$Ca/Bo$) causes the free surface shape to become out of phase with the undulating boundary and leads to a corresponding increase in 
$Q$.
Similar observations can be made by systematically plotting free surface shapes for varying either 
$Re$ or 
$Ca/Bo$, respectively, while the other parameter is fixed. As shown in figure 9(a), we observe that the free surface becomes more deformed and more out of phase with the undulating surface as we vary 
$Re$ from 0 to 1 at 
${Ca}/{Bo}=0.5$. As demonstrated in figure 8(c), this change in the free surface shape corresponds a monotonic increase in 
$Q$, in a manner that is qualitatively similar to the GW case. By contrast, figure 9(b) includes the free surface shapes at 
${Re}=0.4$, as 
$Ca/Bo$ is increased from 0.2 to 1. In this range of parameters, the free surface shape is shown to gradually evolve from flat to more in phase with the undulating surface with a larger amplitude. However, the resultant flux 
$Q$ plotted in figure 8(a) varies non-monotonically, reaching a maximum around 
${Ca}/{Bo}\approx 0.4$ before slightly decreasing with 
$Ca/Bo$. This suggests an inherent shortcoming of determining an ‘optimal’ free surface shape by simply plotting 
$h_2$. To help resolve the unclear connection between the free surface shape and the flux, we will derive the analytical solutions of 
$h_2$ and 
$Q$ in § 4.3.
(a) The theoretical predictions of 
$h_2$ for 
${Ca/Bo}=0.5$ under different 
${Re}$. The different solid lines (coloured from light to dark green) correspond to the values of 
${Re}=0, 0.25, 0.5, 0.75$ and 
$1$. (b) The theoretical predictions of 
$h_2$ for 
${Re}=0.4$ under different 
${Ca/Bo}$. The different solid lines (coloured from light to dark green) correspond to the values of 
${Ca/Bo}=0.2, 0.4, 0.6, 0.8$ and 
$1$. The grey grid lines have been added to (b) to highlight the changes in 
$h_2$ with 
$Ca/Bo$.
4.3. Analytical solutions of $Q$ and $h_2$
To elucidate the mixed effects of 
${Ca}/{Bo}$ and 
$Re$, we obtain an analytical expression for 
$Q$ via asymptotic expansions. First, we recall 
$q=q(h_1,h_2)$ via 
$q =\int _{h_1}^{h_2}\sum _{j=0}^{2} u_j(x,y)\, {\rm d} y$, which contains a term with the prefactor of 
$Bo/Ca$ and that with 
$1/Ca$. Since the magnitude of 
$Bo$ is approximately 
$O(10^2)$ to 
$O(10^3)$ based on the characteristic parameters, we neglect the term prefaced with 
$1/Ca$. Then, we expand 
$q$ as 
$q=q_0+a q_1+a^2 q_2 + a^3q_3+O(a^4)$ by applying an asymptotic expansion of 
$h_2$ in orders of 
$a$
\begin{equation} h_2=1+a f_1+a^2 f_2+ a^3f_3+O(a^4).\end{equation}
In addition, 
$q$ is also a function of the shape of the undulating surface, 
$h_1=\epsilon \sin (2{\rm \pi} x)$, where 
$\epsilon \equiv \delta /H$ and 
$\delta$ is the characteristic undulation amplitude of the solid boundary (Lee et al. Reference Lee, Bush, Hosoi and Lauga2008). Based on the characteristic values from the experiments, we find that 
$a$ and 
$\epsilon$ have the same order of magnitude. Hence, we define an 
$O(1)$ parameter 
$\alpha$, where 
$\epsilon =\alpha a$, and replace 
$\epsilon$ with 
$\alpha a$ to ensure that 
$q$ can be consistently expanded in orders of 
$a$.
The resultant expressions of 
$q$ at each order as functions of the free surface profile are listed in Appendix B, from which we solve for 
$q_0$, 
$q_1$, 
$q_2$ and 
$q_3$, together with 
$f_1$, 
$f_2$ and 
$f_3$. The analytical solutions for 
$q_0$, 
$q_1$, 
$q_2$ and 
$q_3$ and our detailed solution procedure are also included in Appendix B. Once 
$q$ is known at each order, we use 
$Q=q+1$ (i.e. 
$Q_0=1+q_0$, 
$Q_2=q_2$ and 
$Q_3=q_3$) to find the analytical expression of the laboratory-frame flux 
$Q$
\begin{align} Q= \frac{6{\rm \pi}^2 a^2 \alpha^2}{9({Ca/Bo})^2+4{\rm \pi}^2}+\frac{9a^3}{40} \alpha^2\widetilde{Re}\left\{\frac{45({Ca/Bo})^3{\rm \pi}^2 +148({Ca/Bo}){\rm \pi}^4}{(9({Ca/Bo})^2+4{\rm \pi}^2)^2}\right\}+O(a^4),\end{align}
where 
$Q_1=0$. We plot the first two non-zero orders of 
$Q$ in (4.5) as a function of 
$Ca/Bo$ for given 
$Re$, or 
$a \widetilde {Re}$, where 
$a = 0.126$. The resultant plot in figure 10(a) qualitatively matches the full numerical solutions of 
$Q$ in figures 8(a) and 8(b) for the same parameter range. However, the asymptotic results do not exhibit the reduction in 
$Q$ with increasing 
$Re$ in the limit of 
${Ca/Bo}\rightarrow 0$, as this analysis does not include the higher-order viscous terms that give rise to a slight dip in 
$Q$ in this limit. As 
$Re$ and 
$\widetilde {Re}$ simply differ by a constant prefactor, we interpret our asymptotic results in terms of 
$Re$ instead of 
$\widetilde {Re}$ to be consistent with the results shown in figures 7 and 8.
(a) The analytical solutions of 
$Q$ plotted against 
${Ca}/{Bo}$ with different values of 
${Re}$ from 0 to 0.8. (b) Same asymptotic solutions with the range of 
${Ca}/{Bo}$ from 0 to 1.5, with the grey grid lines. The red dashed line depicts the maximum 
$Q$ for any given 
$Re$.
The first term in (4.5) corresponds to 
$a^2Q_2$ and only consists of 
$Ca/Bo$ and monotonically decreases with increasing 
$Ca/Bo$, while the second term, 
$a^3Q_3$, is the product of 
$Re$ and a rational function of 
$Ca/Bo$. Specifically, for given 
$Ca/Bo$, the second term increases linearly with 
$Re$, which validates the effects of inertia to enhance the fluid flux. On the other hand, the effects of 
$Ca/Bo$ on 
$Q$ strongly depend on the value of 
$Re$, which helps determine the relative importance of the two terms. For instance, at 
${Re}=0$, only the first term is retained and leads to a decrease of 
$Q$ with 
$Ca/Bo$, matching the results of the SO data.
To further elucidate the coupled effects of 
$Re$ on 
$Ca/Bo$, we set 
$\partial Q/ \partial ({Ca/Bo})=0$ to compute the critical value of 
$Ca/Bo$ (i.e. 
$Ca/Bo_{cr}$) that yields the maximum value of 
$Q$ at given 
$Re$. The resultant 
$Ca/Bo_{cr}$ is plotted as a red dashed line in figure 10(b). The plots shows that 
$Ca/Bo_{cr}$ monotonically increases with 
$Re$, such that 
$Ca/Bo_{cr}\approx 0.7$ around 
${Re}\approx 0.4$. This explains why, in the range of 
$Ca/Bo$ from 0 to 0.7, we observe the transition of 
$Q$ from monotonically decreasing, being non-monotonic, then to monotonically increasing with 
$Ca/Bo$, as 
$Re$ is increased (see figures 8 and 10a).
As shown in figure 10(b), once we expand the range of 
$Ca/Bo$, we can clearly observe that 
$Q$ always varies non-monotonically with 
$Ca/Bo$ for non-zero 
$Re$. Increasing 
$Re$ simply extends the range of 
$Ca/Bo$ over which 
$Q$ increases, coinciding with larger 
$Ca/Bo_{cr}$. Furthermore, in the limit of 
${Ca/Bo}\gg 1$, we find that 
$Q\approx (6/9){\rm \pi} ^2 a^2\alpha ^2 {(Ca/Bo)^{-2}}+(1/8){\rm \pi} ^2 a^3\alpha ^2 \widetilde {Re}{(Ca/Bo)^{-1}}$, which shows a monotonic decrease with 
$Ca/Bo$ for all 
$Re$. Therefore, we reasonably conclude that 
$Q$ is reduced in the regime where 
$Ca/Bo$ dominates regardless of 
$Re$, while in the intermediate regime where 
$Re$ and 
$Ca/Bo$ are comparable, 
$Q$ is shown to increase with 
$Ca/Bo$.
4.3.1. Connecting free surface shapes to $Q$
Going back to the asymptotic expansions in § 4.3, the free surface shape is given by 
$h_2=1+a f_1+a^2 f_2 + O(a^3)$ in response to the undulating surface of 
$h_1(x)=a\alpha \sin {(2{\rm \pi} x)}$. In the limit of 
$a\ll 1$, we obtain the analytical expressions for 
$f_1$ and 
$f_2$ from (B3) and (B4) in Appendix B
\begin{equation} f_1 = \frac{6{\rm \pi}\alpha({Ca}/{Bo})\cos{(2{\rm \pi} x)} + 9\alpha({Ca}/{Bo})^2\sin{(2{\rm \pi} x)}}{9({Ca}/{Bo})^2 + (2{\rm \pi})^2}, \end{equation}
\begin{align} f_2 &= 12{\rm \pi}^2({Ca/Bo})\alpha\widetilde{Re}\left[\frac{12{\rm \pi}({Ca/Bo})\cos{(2{\rm \pi} x)} + (9({Ca/Bo})^2-4{\rm \pi}^2)\sin{(2{\rm \pi} x)}}{5(9({Ca/Bo})^2+4{\rm \pi}^2)^2}\right]\nonumber\\ &\quad + Q_2({Ca/Bo}) + F(x,{Ca/Bo}), \end{align}
where 
$F(x,{Ca/Bo})$ corresponds to
\begin{align} F &= 6{\rm \pi}^2\alpha^2
\frac{64{\rm \pi}^4 + 180{\rm \pi}^2({Ca/Bo})^2 +
81({Ca/Bo})^4}{(9({Ca/Bo})^2+4{\rm \pi}^2)^2({Ca/Bo})^2+16{\rm \pi}^2)}
+\frac{6{\rm \pi}^2\alpha^2}{Bo/Ca}\nonumber\\ &\quad \times \left[\frac{9({Ca/Bo})(9({Ca/Bo})^2
-20{\rm \pi}^2)\cos{(4{\rm \pi} x)}-24{\rm \pi}(9({Ca/Bo})^2 -
2{\rm \pi}^2)\sin{(4{\rm \pi}
x)}}{(9({Ca/Bo})^2+4{\rm \pi}^2)^2({Ca/Bo})^2+16{\rm \pi}^2)}\right].
\end{align}
Notably, 
$f_1$ is a linear superposition of 
$\cos {(2{\rm \pi} x)}$ and 
$\sin {(2{\rm \pi} x)}$, independent of 
$Re$, while 
$f_2$ comprises sine and cosine functions of higher wavenumbers and exhibits dependence on both 
$Re$ and 
$Ca/Bo$. Specifically, the periodic functions with higher wavenumbers in 
$f_2$ suggest the effects of nonlinearities in our model. However, since they persist even at 
$\widetilde {Re}=0$, they must stem from the nonlinearity of the free surface itself and are not related to the inertia term in the governing equation (i.e. (3.6)). On the other hand, fluid inertia from (3.6) leads to the term in 
$f_2$ that is linear in 
$\widetilde {Re}$. As will be further discussed, the 
$\widetilde {Re}$ term causes the free surface to become out of phase with 
$h_1$, which helps enhance the fluid flux.
Figure 11 includes the plots 
$f_1$ and 
$f_2$ for different dimensionless parameters; note that the dashed curve corresponds to the shape of the undulating surface, or 
$h_1/a = \alpha \sin {(2{\rm \pi} x)}$, as a reference point. As shown in figure 11(a), 
$f_1$ is uniformly zero at 
${Ca}/{Bo}=0$ and gradually deforms more with increasing 
$Ca/Bo$. Then, as 
${Ca}/{Bo}$ reaches 10 (see the inset of figure 11a), 
$f_1$ closely matches the shape of the undulating surface, as the 
$\alpha \sin {(2{\rm \pi} x)}$ term comes to dominates in the limit of 
${Ca/Bo}\rightarrow \infty$.
(a) The profiles of 
$f_1$ under different 
$Ca/Bo$. The different solid lines (coloured from light to dark green) correspond to the different values of 
${Ca/Bo}=0, 0.25, 0.5, 0.75$ and 
$1$. The inset shows the profiles of 
$f_1$ for 
${Ca/Bo}=1-10$. (b) The profile of 
$f_2$ at 
${Re}=0.5$ for 
${Ca/Bo}=0, 0.5, 1, 1.5$ and 
$2$ and (c) for 
${Ca/Bo}=2, 2.5, 3, 3.5$ and 
$4$. (d) The profile of 
$f_2$ at 
${Ca/Bo}=0.5$ under different 
${Re}$. The different solid lines (coloured from light to dark green) correspond to the different values of 
${Re}=0, 0.25, 0.5, 0.75$ and 
$1$. The red-dashed lines in the figures correspond to the shape of 
$\alpha \sin (2{\rm \pi} x)$.
For 
$f_2$ that depends on both 
$Re$ and 
$Ca/Bo$, we first consider the effects of 
$Ca/Bo$ on 
$f_2$ while 
$Re$ is held constant at 0.5. Similar to the evolution of 
$f_1$, 
$f_2$ becomes more deformed as 
$Ca/Bo$ is increased from 0 to 2, as shown in figure 11(b). However, distinct from 
$f_1$, 
$f_2$ has a flattened region near 
$x=0$ that grows in amplitude with 
$Ca/Bo$, due to the contributions from constants and higher wavenumber terms in (4.8). Then, as 
$Ca/Bo$ is increased above 2 (see figure 11c), the flat region near 
$x=0$ appears to saturate and no longer grow with 
$Ca/Bo$, while the deformation of 
$f_2$ away from 
$x=0$ decreases. This non-monotonic behaviour of 
$f_2$ with 
$Ca/Bo$ may have connections to the non-monotonicity in the flux 
$Q$, which will be further explored. Finally, when we instead fix 
$Ca/Bo$ at 0.5 and vary 
$Re$, figure 11(d) shows that the wavenumber of 
$f_2$ is reduced with increasing 
$Re$, while the amplitude of deformation increases. Furthermore, 
$f_2$ becomes more out of phase with the undulating surface for increasing 
$Re$, which suggests the formation of a slight dip and a corresponding increase in 
$Q$ as observed in the GW mixture case.
To more clearly understand the connection between the free surface shapes and 
$Q$, we write the resultant flux explicitly as a function of 
$f_1$ and 
$f_2$, such that
\begin{gather} Q_2 = \frac{Bo}{Ca}\int^1_0f_1'(f_1 -\alpha \sin{(2{\rm \pi} x)})\, {\rm d}\kern0.06em x, \end{gather}
\begin{gather}Q_3= \frac{Bo}{Ca}\int^1_0[f_2'(f_1 - \alpha\sin{(2{\rm \pi} x)}) + \widetilde{Re}f_1'(16{\rm \pi}\alpha\cos{(2{\rm \pi} x)} - 5f_1')/48 + f_1'f_2]\, {\rm d}\kern0.06em x, \end{gather}
which have been reduced from (B3) and (B4) via integration (i.e. 
$\int _{0}^{1} Q_2 \, {\rm d}\kern0.06em x=Q_2$ and 
$\int _{0}^{1} Q_3 \, {\rm d}\kern0.06em x=Q_3$). Here, the prime denotes the differentiation with respect to 
$x$. Note that combining (4.7)–(4.10) leads to the expression for 
$Q$ given in (4.5). Clearly, 
$Q_2$ scales with 
$f_1 - \alpha \sin {(2{\rm \pi} x)}$, which explains the reduction in 
$Q$ as 
$f_1$ approaches the shape of the undulating surface with increasing 
$Ca/Bo$. Furthermore, 
$Q_2$ is completely independent of 
$Re$ and shows a monotonic decrease with 
$Ca/Bo$.
The effects of inertia on the flux are introduced exclusively through 
$Q_3$ in two ways: via terms that are linear in 
$Re$ and those that depend on 
$f_2$. Hence, we divide up 
$Q_3$ in the following way:
\begin{gather} Q_{30}= \frac{\widetilde{Re}}{48Ca/Bo}\int^1_0f_1' (16{\rm \pi}\alpha\cos{(2{\rm \pi} x)} - 5f_1')\, {\rm d}\kern0.06em x = \frac{9{\rm \pi}({Ca/Bo})\alpha^2\widetilde{Re}}{72 ({Ca/Bo})^2+32{\rm \pi}^2}, \end{gather}
\begin{gather}Q_{31} = \frac{Bo}{Ca}\int^1_0[f_2'(f_1 - \alpha\sin{(2{\rm \pi} x)}) + f_1'f_2]\, {\rm d}\kern0.06em x = \frac{144{\rm \pi}^4({Ca/Bo})\alpha^2\widetilde{Re}}{5(9 ({Ca/Bo})^2+4{\rm \pi}^2)^2}, \end{gather}
so that 
$Q_{30}$ contains the effects of inertia from the changes in the flow itself, while 
$Q_{31}$ includes the effects of inertia through the changes in the free surface shape via 
$f_2$. Interestingly, the higher wavenumber terms in 
$f_2$ drop out of 
$Q_{31}$ upon integration, so that only inertia-induced changes in the free surface contribute to the fluid flux at 
$O(a^3)$.
We plot 
$Q_{30}$ (blue) and 
$Q_{31}$ (red) together in figure 12 as a function of 
$Ca/Bo$ for 
${Re} = 0.5$ and also for varying 
$Re$ for 
${Ca}/{Bo} = 0.5$, respectively. The results in figure 12 demonstrate that 
$Q_{31}$ is overall larger than 
$Q_{30}$. Hence, the effects of inertia enter more strongly through the changes in the free surface shape than through the inertia-induced changes in the flow. Specifically, figure 12(b) shows a monotonic increase in both 
$Q_{30}$ and 
$Q_{31}$ for increasing 
$Re$. This increase in 
$Q_{31}$ can be understood via the gradual increase in 
$f_2$ that becomes more out of phase with 
$h_1$ as 
$Re$ is increased (see figure 11d).
(a) Values of 
$Q_{30}$ and 
$Q_{31}$ plotted against 
$Ca/Bo$ with 
${Re}=0.5$ and (b) 
$Q_{30}$ and 
$Q_{31}$ plotted against 
$Re$ with 
${Ca/Bo}=0.5$.
By contrast, 
$Q_{30}$ and 
$Q_{31}$ in figure 12(a) show non-monotonic dependence on 
$Ca/Bo$ at 
${Re}=0.5$, reaching a maximum at 
${Ca}/{Bo} = 2{\rm \pi} /3$ and 
$2{\rm \pi} /(3\sqrt {3})$, respectively. Curiously, the transitional behaviour in 
$f_2$ as observed in figures 11(b) and 11(c) occurs around 
${Ca}/{Bo} \approx 2$, which is larger than the local maximum of 
$Q_{31}$ (i.e. 
${Ca}/{Bo} = 2{\rm \pi} /(3\sqrt {3})$). This difference in the transitional points for 
$f_2$ and 
$Q_{31}$ may be explained by the fact that 
$Q_{31}$ also contains a term that scales with 
$f_1-\alpha \sin {(2{\rm \pi} x)}$. As 
$f_1-\alpha \sin {(2{\rm \pi} x)}$ monotonically decreases towards zero with 
$Ca/Bo$, this term mitigates the effects of 
$f_2$ on the flux and causes 
$Q_{31}$ to reach a maximum at lower 
$Ca/Bo$. Overall, we find that the changes in the flux with 
$Re$ and 
$Ca/Bo$ can be sufficiently explained through the corresponding changes in the free surface shape.
4.3.2. Selection of the working fluid
While decoupling 
$Re$ and 
$Ca/Bo$ is crucial in isolating the effects of inertia, in practice, we select a working fluid (i.e. 
$\rho$, 
$\mu$, 
$\sigma$) with given film thickness 
$H$ and then systematically vary 
$V_w$, which causes 
$Re$ and 
$Ca/Bo$ to change simultaneously. As we have experimentally demonstrated in § 4.1, once the working fluid has been selected, the resultant pump rate 
$Q$ may increase or decrease with 
$V_w$, which we have rationalised in terms of 
$Re$ and 
$Ca/Bo$ and with corresponding free surface shapes. Building on our original experimental observations, we will discuss how the selection of working fluids determines the dependence of 
$Q$ on 
$V_w$.
Going back to the phase map in figure 7, we first seek a boundary that separates 
$Q$ that either increases or decreases with 
$V_w$, where the boundary satisfies the condition: 
${\rm d}Q/{\rm d}V_w=0$. In addition, we recall that 
$Q$ for a given fluid can be represented by a line passing through the origin of the phase map, whose slope is given by 
${Re}/({Ca}/{Bo})\propto (\rho /\mu )^2$ from (4.3). For instance, on the phase diagram in figure 7, SO corresponds to a straight line with a slope of 
${Re}/({Ca}/{Bo})=0.039$, while the GW mixture has a slope of 2.983. Hence, we will analyse how the boundary of 
${\rm d}Q/{\rm d}V_w=0$ interacts with straight lines that represent unique fluids. As the analytical solution of 
$Q$ in (4.5) is a function of 
$Re$ and 
$Ca/Bo$, we apply the chain rule
\begin{equation} \frac{{\rm d}Q}{{\rm d}V_w}=\frac{\partial Q}{\partial Re} \frac{\partial Re}{\partial V_w}+\frac{\partial Q}{\partial (Ca/Bo)} \frac{\partial (Ca/Bo)}{\partial V_w}=0,\end{equation}so that
\begin{align} &\frac{9a^3 \widetilde{Re}\alpha^2}{40}\left\{\frac{45({Ca/Bo})^2{\rm \pi}^2 +148{\rm \pi}^4}{(9({Ca/Bo})^2+4{\rm \pi}^2)^2}\right\}+\frac{9a^3 \alpha^2 \widetilde{Re}(135({Ca/Bo})^2{\rm \pi}^2 +148{\rm \pi}^4)}{40(9({Ca/Bo})^2+4{\rm \pi}^2)^2}\nonumber\\ &\qquad -\frac{108a^2 \alpha^2 ({Ca/Bo}){\rm \pi}^4}{(9({Ca/Bo})^2+4{\rm \pi}^2)^2}-\frac{81a^3 \alpha^2 \widetilde{Re}(45({Ca/Bo})^4{\rm \pi}^2 +148({Ca/Bo})^2{\rm \pi}^4)}{10(9({Ca/Bo})^2+4{\rm \pi}^2)^3}=0\nonumber\\ &\quad \equiv G(\widetilde{Re},Ca/Bo). \end{align} We plot 
$G(\widetilde {Re},Ca/Bo)$ as a dashed line in figure 13(a), which corresponds to the boundary of 
${\rm d}Q/{\rm d}V_w=0$. The boundary has a positive slope that monotonically increases with 
$Ca/Bo$, so the minimum slope at the origin is given by 
${Re}/({Ca}/{Bo})\approx 0.164$. Therefore, a fluid with a slope smaller than 0.164 will yield 
$Q$ that monotonically decreases with 
$V_w$. On the other hand, the critical slope that separates non-monotonic 
$Q$ from increasing 
$Q$ depends on the range of 
$V_w$ considered. For instance, the critical value of 
$Re/(Ca/Bo)$ above which 
$Q$ increases monotonically with 
$Ca/Bo$ corresponds to 0.243, in the current range of 
$Re$ from 0 to 1. Hence, for a fluid with 
$Re/(Ca/Bo)$ larger than 0.164, the behaviour of the resultant 
$Q$ with 
$V_w$ depends on by how much the slope exceeds 0.164 and the range of 
$V_w$ considered.
(a) The phase diagram summarising the magnitude of 
${\rm d}Q/{\rm d}V_w$ with varying 
$Re$ and 
$Ca/Bo$. The black dot-dashed line refers to 
${\rm d}Q/{\rm d}V_w=0$, with the grey contour lines corresponding to constant 
${\rm d}Q/{\rm d}V_w$. (b) Three lines whose slopes are given by 
$Re/(Ca/Bo)= 0.141$ (blue), 0.188 (cyan) and 0.243 (green) exhibit distinct behaviours of 
$Q$ with 
$V_w$.
To demonstrate this, in figure 13(b), we have plotted 
$Q$ as a function of 
$V_w$ for three different values of 
$Re/(Ca/Bo)$, representing three different fluid properties. Specifically, for a fluid whose slope is slightly higher than 0.164 (i.e. 
$Re/(Ca/Bo) = 0.188$), 
$Q$ first increases with 
$V_w$ and then begins to decrease with 
$V_w$ once the line intersects 
$G(\widetilde {Re},Ca/Bo)$, shown as cyan in figure 13(b). When 
$Re/(Ca/Bo)$ is set even higher at 0.243 (green line), 
$Q$ monotonically increases with 
$V_w$, in a qualitatively similar manner as the GW mixture. Finally, for the fluid with 
$Re/(Ca/Bo) = 0.141 < 0.164$ shown in a solid blue line, 
$Q$ monotonically decreases with 
$V_w$, consistent with the result of SO. Hence, we have developed a simple criterion for predicting how 
$Q$ will vary with 
$V_w$ based on the fluid properties and the film geometry.
5. Conclusion
In summary, we experimentally and theoretically analyse the dynamics of thin, free surface flows that are driven by an undulating surface. We experimentally demonstrate that the net pumping rate 
$Q$ of the free surface flow behaves in a qualitatively different way depending on 
$Re$, which cannot be explained by the standard lubrication model of Pandey et al. (Reference Pandey, Chen, Yuk, Sun, Roh, Takagi, Lee and Jung2023). To rationalise the experimental observations, we construct a two-dimensional mathematical model that incorporates the effects of finite inertia by applying an asymptotic expansion on the velocity and pressure fields. Our theoretical model correctly captures the behaviour of 
$Q$ for varying wave speed 
$V_w$ with no fitting parameters. The theoretical model also provides the prediction of the free surface profile for increasing 
$V_w$. Specifically, in the viscosity-dominated regime, we find that the free surface profile gradually conforms to the shape of the undulating surface and yields a decrease in 
$Q$, as observed in the experiments with SO. By contrast, in the regime with finite inertia, our model predicts that the free surface becomes out of phase with the bottom surface, which leads to a corresponding increase in 
$Q$. This is consistent with the experiments of the GW mixture.
In addition to rationalising experimental observations, we construct a phase map with numerical simulations to decouple the effects of inertia 
$(Re)$ from viscous and hydrostatic effects 
$(Ca/Bo)$. The simulation results show a decrease in 
$Q$ with 
$Ca/Bo$ for small 
$Re$, but as 
$Re$ is increased, 
$Q$ transitions to increasing with 
$Ca/Bo$, suggesting the potential mixed effects of 
$Ca/Bo$ and 
$Re$. To further delve into these mixed effects, we derive the analytical expressions of 
$Q$ and the free surface shape 
$h_2$, by expanding them in the limit of asymptotically small deformations. The resultant analytical expression of 
$Q$ yields a critical value of 
$Ca/Bo_{cr}$ above which 
$Q$ decreases with 
$Ca/Bo$. Specifically, 
$Ca/Bo_{cr}$ increases with 
$Re$, demonstrating that inertia tends to enhance 
$Q$ and delays the dominance of the viscous effects over hydrostatics that reduce 
$Q$. We also obtain the connection between the free surface shapes 
$h_2$ and 
$Q$ by analysing the analytical expressions of 
$h_2$. Finally, by setting 
${\rm d}Q/{\rm d}V_w=0$, we demonstrate how the choice of the working fluid directly influences the dependence of 
$Q$ on 
$V_w$.
Overall, our two-dimensional mathematical model successfully captures the dynamics of free surface flows with an undulating surface. However, the model incorporates a number of major assumptions that could be improved upon in the future studies. For example, as 
$Re$ is further increased, our preliminary experiments have demonstrated the generation of free surface waves that are not dictated by that of the undulating surface. Understanding this new flow regime will require bringing in the time dependence into our model, which we plan to develop in the future. In addition, due to the size of the undulator and corresponding large 
$Bo$, our current analysis shows negligible surface tension effects in determining the pumping rate. We hypothesise that surface tension will come to play a dominant role in the limit of small 
$Bo$, which we are interested in exploring with a miniaturised undulator. Finally, our future experimental studies may include adding particles on the undulator-driven free surface and exploring the effects of inertia on their trajectories (Leal Reference Leal1980).
Funding
This work was partly funded by the National Science Foundation (CMMI-2042740).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical scheme
We first set an initial guess for the flux 
$q=0.8$ (step (i)) and start the simulation with 
$V_w=0.01 \ {\rm m}\ {\rm s}^{-1}$. Since the free surface deformation should be minor for small 
$V_w$, we set a flat free surface (
$h_2=1$) as the initial guess for 
$h_2$ and use the MATLAB routine ‘bvp4c’ to compute 
$h_2$, subject to the periodic boundary conditions (step (ii)). Next, we check whether or not the total fluid area 
$A_i$ for the converged solution of 
$h_2$ is conserved, or 
$A_i = 1$ (step (iii)). If 
$A_i\ne 1$, we incrementally change the value of 
$q$ and repeat the process until 
$A_i=1$ (step (iv)). Note that we reduce 
$q$ if 
$A_i>1$, whereas we increase 
$q$ if 
$A_i<1$. In addition, to ensure that the simulation is stable, we apply the continuation method and use the previous solution of 
$h_2$ as the initial guess, as we systematically increase 
$V_w$.
For generating the solutions as a function of 
$Ca/Bo$ and 
$Re$ (§ 4.2), we follow the same numerical scheme as described above. Specifically, to simultaneously obtain the flux and free surface shape for varying 
$Ca/Bo$, we first set the value of 
$Re$ and start the simulation with 
${Ca}/{Bo}=0.01$ with a flat free surface (
$h_2=1$) as the initial guess. Then, following the same steps as above, we obtain 
$h_2$ and 
$q$ for systematically increasing 
$Ca/Bo$. For varying 
$Re$, we first fix the value of 
$Ca/Bo$ and start the simulation with 
${Re}=0.01$. However, the initial profile of 
$h_2$ is not set to be flat for large 
$Ca/Bo$. Instead, we replace 
$h_2=1$ with 
$(h_2, h_2', h_2'', h_2''', h_2'''', h_2''''')=(1, 0.1, 0, 0, 0)$ and change the last three guesses with small adjustments (e.g. 0.01), to obtain the converged solution of 
$h_2$ via ‘bvp4c’ if necessary.
Appendix B. Asymptotic expansions of $q$
Expanding 
$q$ via the asymptotic expansion of 
$h_2$ yields the following expressions for 
$q_0$, 
$q_1$, 
$q_2$ and 
$q_3$:
\begin{gather} q_0={-}1, \end{gather}
\begin{gather}q_1={-}f_1+ \alpha \sin(2{\rm \pi} x)-\frac{1}{3 {Ca/Bo}}f'_1, \end{gather}
\begin{gather}q_2={-}f_2 - \frac{1}{Ca/Bo} f_1 f_1' + \frac{\alpha}{Ca/Bo} \sin(2{\rm \pi} x) f_1'+\frac{1}{3Ca/Bo}f_2'-\frac{\widetilde{Re}}{15 Ca/Bo}f_1'', \end{gather}
\begin{align} q_3&={-}f_3+\frac{{\rm \pi}\alpha \widetilde{Re}}{3 Ca/Bo} \cos(2{\rm \pi} x) f_1'-\frac{1}{Ca/Bo}f_1'(f_1^2-f_2)+\frac{2\alpha}{Ca/Bo}\sin(2{\rm \pi} x)f_1f_1''\nonumber\\ &\quad -\frac{\alpha^2}{Ca/Bo}\sin^2(2{\rm \pi} x)f_1 \frac{5\widetilde{Re}}{48Ca/Bo}f_1'^2-\frac{1}{Ca/Bo}f_1 f_2'+ \frac{\alpha}{Ca/Bo}\sin(2{\rm \pi} x)f_2'-\frac{1}{3Ca/Bo}f_3'\nonumber\\ &\quad-\frac{\widetilde{Re}}{3Ca/Bo}f_1 f_1''+ \frac{\widetilde{Re}\alpha}{3Ca/Bo}\sin(2{\rm \pi} x)f_1''- \frac{37\widetilde{Re}\alpha}{1680(Ca/Bo)^2}f_1' f_1''- \frac{\widetilde{Re}}{15Ca/Bo} f_2''\nonumber\\ &\quad-\frac{1}{40Ca/Bo} f_1'''-\frac{17\widetilde{Re}^2}{1260Ca/Bo} f_1'''. \end{align}
We note that 
$q_0=-1$ as a consequence of being in the wave frame.
Next, we simultaneously solve for 
$q$ and 
$f$ at each order via the ‘DSolve’ function in Mathematica, by imposing the periodic boundary conditions and the conservation of fluid area. Note that 
$\int _{0}^{1} (h_2-h_1)\, {\rm d}\kern0.06em x=1$ can be simplified to 
$\int _{0}^{1} h_2 \, {\rm d}\kern0.06em x=1$, since 
$h_1=\alpha a\sin (2{\rm \pi} x)$. Also, the first term in the expansion of 
$h_2$ is 1, which automatically satisfies the total fluid area. Hence, the conservation of fluid area for the remaining terms in the expansion becomes 
$\int _{0}^{1} f_n \, {\rm d}\kern0.06em x=0$ for 
$n=1,2,3$. The resultant analytical solutions for 
$q_1$, 
$q_2$ and 
$q_3$ correspond to
\begin{gather} q_1=0, \end{gather}
\begin{gather}q_2=\frac{6{\rm \pi}^2\alpha^2}{9({Ca/Bo})^2+4{\rm \pi}^2}, \end{gather}
\begin{gather}q_3=\frac{9\widetilde{Re}(45{\rm \pi}^2({Ca/Bo})^3+148 {\rm \pi}^4({Ca/Bo}))}{40(9({Ca/Bo})^2+4{\rm \pi}^2)^2}. \end{gather}
Open access
