2.1. Problem formulation

The problem of interest is that of a gravity-driven layer of Newtonian fluid, with constant density $\rho$ and dynamic viscosity $\mu$, flowing over a uniformly heated, periodically corrugated rigid substrate inclined at an angle $\beta$ to the horizontal, as illustrated schematically in figure 1. The corresponding Cartesian coordinate system is orientated such that the $X$-axis points along and down the inclined substrate which is considered to be infinite and invariant in the $Y$-direction, with the $Z$-axis normal to it, rendering the principal problem two-dimensional. The substrate profile, $S(X)$, is measured relative to the $X$-axis and given by

(2.1)\begin{equation} S(X) = A \frac{ \left[ 1 + \cos \left( \dfrac{ 2 {\rm \pi}X }{ L } \right) \right] }{2}, \end{equation}

where $A$ is the corrugation amplitude and $L$ its wavelength.

Figure 1.

Schematic of film flow down a uniformly heated, wavy rigid substrate inclined at an angle $\beta$ to the horizontal, showing the main geometrical features relative to the adopted coordinate system.

The film thickness, $H(X,T)$, at a downstream location $X$, is the difference between the free-surface location, $Z = F(X,T)$, and the corrugation height, $Z = S(X)$. The temperature of the substrate, $\varTheta _s$, and that of the surrounding ambient gas, $\varTheta _a$, remain fixed and constant with the difference between them defined as $\varTheta _{\varDelta } = \varTheta _s - \varTheta _a$.

The equations governing the flow are the continuity and Navier–Stokes (NS) equations

(2.2a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}$$

(2.2b)$$\begin{gather}\rho \frac{ \textrm{D} \boldsymbol{U} }{ \textrm{D} T } ={-} \boldsymbol{\nabla} P + \rho \boldsymbol{{g}} + \mu \nabla^2 \boldsymbol{U}, \end{gather}$$

where $\boldsymbol {U} = (U , W)^{\textrm {T}}$ and $P$ are the velocity and pressure fields, respectively, $\boldsymbol {\nabla } = ( \partial / \partial X , \partial / \partial Z )$ is the differential operator and $\boldsymbol {g} = \bar {g} ( \sin {\beta } , - \cos {\beta })^{\textrm {T}}$ is the acceleration due to gravity (with magnitude $\bar {g}$); together with a corresponding convection–diffusion equation representing the conservation of energy. For the case of no volumetric heating and negligible thermal viscous dissipation, the latter can be written in terms of the temperature field $\varTheta$ as

(2.3)\begin{equation} \rho C_P \frac{D \varTheta}{D T} = \kappa \nabla^2 \varTheta, \end{equation}

in which the heat capacity at constant pressure $C_P$ and thermal conductivity $\kappa$ of the fluid are both assumed to remain constant.

The boundary conditions at the non-porous substrate include a no-slip and constant temperature requirement, namely

(2.4a,b)\begin{equation} \boldsymbol{U} = 0,\quad \varTheta = \varTheta_s, \quad \text{at} \ Z = S(X). \end{equation}

The fluid is assumed to be non-volatile and $\varTheta _{\varDelta }$ to be small enough that no evaporation of the liquid film takes place. Accordingly, the evolution of the free surface is described by the kinematic condition, given by

(2.5)\begin{equation} \frac{ \partial H }{ \partial T } + U \frac{ \partial F }{ \partial X } = W, \quad \text{at} \ Z = F(X,T), \end{equation}

where $F(X,T) = S(X) + H(X,T)$. The normal and tangential force balance equations there are

(2.6a,b)\begin{equation} \left. \begin{array}{ll} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \hat{\boldsymbol{n}} - \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{\hat{T}} \boldsymbol{\cdot} \hat{\boldsymbol{n}} ={-} \sigma ( \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{n}})\\ \ \ \hat{\boldsymbol{t}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \hat{\boldsymbol{n}} - \hat{\boldsymbol{t}} \boldsymbol{\cdot} \boldsymbol{\hat{T}} \boldsymbol{\cdot} \hat{\boldsymbol{n}} = \boldsymbol{\nabla} \sigma \boldsymbol{\cdot} \hat{\boldsymbol{t}} \end{array} \right\}\quad \text{at} \ Z = F(X,T), \end{equation}

where $\boldsymbol {T} = - P \boldsymbol {I} + \mu [ ( \boldsymbol {\nabla } \boldsymbol {U} ) + ( \boldsymbol {\nabla } \boldsymbol {U} )^{\textrm {T}} ]$ and $\boldsymbol {\hat {T}} = - P_0 \boldsymbol {I}$ are the stress tensors of the fluid and ambient surrounding air, respectively; $P_0$ is the atmospheric pressure and $\boldsymbol {I}$ is the identity matrix. The normal and tangent unit vectors at the free surface are given by

(2.7a,b)\begin{equation} \hat{\boldsymbol{n}} = \frac{1}{ \sqrt{ 1 + G } } \left( - \frac{ \partial F }{ \partial X } , 1 \right), \quad \hat{\boldsymbol{t}} = \frac{1}{ \sqrt{ 1 + G } } \left( 1 , \frac{ \partial F }{ \partial X } \right); \end{equation}

both of which include the surface curvature pre-factor, $G = ( \partial F / \partial X )^2$. The surface tension of the film, $\sigma$, varies with $\varTheta$ in the following manner:

(2.8)\begin{equation} \sigma = \sigma_0 \left[ 1 + ( \varTheta - \varTheta_a ) \frac{ \partial \sigma }{ \partial \varTheta } \right] , \end{equation}

with $\sigma _0$ the surface tension at temperature $\varTheta _a$; the surface tension gradient $\partial \sigma / \partial \varTheta$ is taken to be constant.

The heat transferred at the free surface follows Newton's law of cooling, which reads

(2.9)\begin{equation} - \kappa \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{\nabla} \varTheta = \alpha (\varTheta - \varTheta_a), \quad \text{at} \ Z=F(X,T), \end{equation}

where $\alpha$ is the liquid–gas heat transfer coefficient there and is valid as long as convective heat transfer is minimal, which can be expected when $\varTheta _{\varDelta }$ is small.

The following scalings are used to non-dimensionalise the governing equations and boundary conditions, (2.1)–(2.9):

(2.10)\begin{equation} \left. \begin{array}{ccc@{}} X = L x, & Z = H_0 z, & T = \dfrac{L}{U_0} t,\\ U = U_0 u, & W = \dfrac{ H_0 U_0 }{ L } w, & \varTheta = \varTheta_{\varDelta} \theta + \varTheta_a,\\ S = H_0 s & H = H_0 h, & P = \dfrac{ \mu U_0 L }{ H_0^2 } p, \end{array} \right\}, \end{equation}

with lowercase letters representing dimensionless variables. The vertical length scale, $H_0$, corresponds to the Nusselt film thickness

(2.11)\begin{equation} H_0 = \left( \frac{ 3 \mu \dot{V_0} }{ \rho \bar{g} \sin{\beta} } \right)^{1/3}, \end{equation}

while the free-surface velocity for a flat film, $U_0 = 3 \dot {V_0} / 2 H_0$, is taken as the velocity scale. The two are linked through the streamwise flow rate per unit cross-sectional width, $\dot {V_0}$.

Applying the scalings leads to the following dimensionless substrate profile:

(2.12)\begin{equation} s(x) = \frac{A}{H_0} \frac{[ 1 + \cos ( 2 {\rm \pi}x )] }{2} \end{equation}

and dimensionless governing equations

(2.13a)$$\begin{gather} \frac{ \partial u }{ \partial x } + \frac{ \partial w }{ \partial z } = 0, \end{gather}$$

(2.13b)$$\begin{gather}{\epsilon} Re \left[ \frac{ \partial u }{ \partial t } + u \frac{ \partial u }{ \partial x } + w \frac{ \partial u }{ \partial z } \right] ={-} \frac{ \partial p }{ \partial x } + 2 + {\epsilon}^2 \frac{ \partial^2 u }{ \partial x^2 } + \frac{ \partial^2 u }{ \partial z^2 },\end{gather}$$

(2.13c)$$\begin{gather}{\epsilon}^3 Re \left[ \frac{ \partial w }{ \partial t } + u \frac{ \partial w }{ \partial x } + w \frac{ \partial w }{ \partial z } \right] ={-} \frac{ \partial p }{ \partial z } - 2 {\epsilon} \cot \beta + {\epsilon}^2 \frac{ \partial }{ \partial x } \left[ {\epsilon}^2 \frac{ \partial w }{\partial x } + \frac{ \partial u }{ \partial z } \right]+ 2 {\epsilon}^2 \frac{ \partial^2 w }{\partial z^2 },\end{gather}$$

(2.13d)$$\begin{gather}{\epsilon} Re\,Pr \left[ \frac{ \partial \theta }{ \partial t } + u \frac{ \partial \theta }{ \partial x } + w \frac{ \partial \theta }{ \partial z } \right] = {\epsilon}^2 \frac{ \partial^2 \theta }{ \partial x^2 } + \frac{ \partial^2 \theta }{ \partial z^2 },\end{gather}$$

featuring the following dimensionless groups: the shallowness parameter, ${\epsilon } = H_0 / L$; the Reynolds number, $Re = \rho U_0 H_0 / \mu$; the Prandtl number, $Pr = ( \mu C_P ) / \kappa$.

The dimensionless boundary conditions at the substrate, $z = s(x)$, read

(2.14a,b)\begin{equation} u|_{z=s} = w|_{z=s} = 0,\quad \theta |_{z=s}= 1, \end{equation}

while scaling and expanding the boundary conditions at the free surface, $z = f(x,t)$, leads to the following expressions:

(2.15a)$$\begin{gather} w|_{z=f} = \frac{ \partial h }{ \partial t } + u|_{z=f} \frac{ \partial f }{ \partial x },\end{gather}$$

(2.15b)$$\begin{gather}p|_{z=f} = p_0 - \left. 2 {\epsilon}^2 \left[ \frac{( 1 - {\epsilon}^2g) \dfrac{ \partial u }{ \partial x } + \dfrac{ \partial f }{ \partial x } \left[ {\epsilon}^2 \dfrac{ \partial w }{ \partial x } + \dfrac{ \partial u }{ \partial z } \right] }{ 1 + {\epsilon}^2 g } \right] \right|_{z=f}- \frac{ {\epsilon}^3 ( 1 - Ma \vartheta ) \dfrac{ \partial^2 f }{ \partial x^2 } }{ Ca [ 1 + {\epsilon}^2g ]^{3/2} }, \end{gather}$$

(2.15c)$$\begin{gather}\left. \left[ {\epsilon}^2 \frac{ \partial w }{ \partial x } + \frac{ \partial u }{ \partial z } \right] \right|_{z=f} = 4 \frac{ \left. {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial u }{ \partial x } \right|_{z=f} }{ 1 - {\epsilon}^2 g } - \left. {\epsilon} \frac{Ma}{Ca} \frac{ \sqrt{ 1 + {\epsilon}^2 g } }{ 1 - {\epsilon}^2 g } \left[ \frac{ \partial \theta }{ \partial x } + \frac{ \partial f }{ \partial x } \frac{ \partial \theta }{ \partial z } \right] \right|_{z=f}, \end{gather}$$

(2.15d)$$\begin{gather}\left. \frac{ \partial \theta }{ \partial z } \right|_{z=f} = \left. {\epsilon}^2 \frac{ \partial f }{ \partial x } \frac{ \partial \theta }{ \partial x } \right|_{z=f} - Bi \vartheta \sqrt{ 1 + {\epsilon}^2 g }. \end{gather}$$

These include the dimensionless free-surface profile, $f(x,t) = s(x) + h(x,t)$, and temperature, $\vartheta = \theta |_{z=f}$, together with the dimensionless surface curvature pre-factor, $g = ( \partial f / \partial x )^2$. In addition, they contain the following additional dimensionless groups: Capillary number, $Ca = \mu U_0 / \sigma _0$; Marangoni number, $Ma = \varTheta _{\varDelta } ( - \partial \sigma / \partial \varTheta )$; Biot number, $Bi = \alpha H_0 / \kappa$.

The fluid pressure is obtained by integrating equation (2.13c) between $z$ and $z = f(x,t)$ with the upper limit of integration given by (2.15b), leading to

(2.16)\begin{align} p &= p_0+ 2 {\epsilon} ( f - z ) \cot \beta - {\epsilon}^2 \frac{ \partial u }{ \partial x } - {\epsilon}^2 \frac{ \partial }{ \partial x } ( u|_{z=f} ) - { \epsilon}^3 \frac{ \partial }{ \partial x } \left[ \frac{ (1 - Ma \vartheta ) }{ Ca } \frac{ \dfrac{ \partial f }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] \nonumber\\ &\quad + {\epsilon}^3 Re \int^{f}_{z} \left[ \frac{ \partial w }{ \partial t } + \frac{ \partial }{ \partial x} ( u w ) + \frac{ \partial }{ \partial z } ( w^2 )\right] \textrm{d}z - {\epsilon}^4 \frac{ \partial }{ \partial x } \int^{f}_{z} \frac{ \partial w }{ \partial x } \, \textrm{d}z, \end{align}

a complete derivation of which is provided in the online supplementary material available at https://doi.org/10.1017/jfm.2021.920.

2.2. Reduced asymptotic model

Solving the full equation set (2.13)–(2.15) represents a daunting task due to the a priori unknown location of the free surface. Accordingly, an asymptotic model of reduced dimensionality is derived following the modelling approach of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000). Proceeding in this way is superior to the gradient expansion of Benney (Reference Benney1966) because it captures the distortion of the primary parabolic flow as deviations manifest.

Gravity-driven film flow is well suited to analysis via a long-wave expansion due to the predominance of surface tension with regard to the film stability, in that long waves are commonly the most unstable modes (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997). Therefore, it is reasonable to assume quantities vary slowly in time and in the $x$-direction, which is to say derivatives with respect to $x$ and $t$ are small. The smallness of space–time derivatives in the present scaling is represented by the shallowness parameter, ${\epsilon } < 1$, with the order of such terms in the long-wave expansion given by the power of the scaling factor, ${\epsilon }^n$. In contrast, gradients with respect to the $z$-coordinate can be considered large. Accordingly, and on the assumption that $u$ and $\theta$ are infinitely differentiable with respect to $z$, the fluid velocity and temperature are expanded as infinite series which take the form

(2.17a,b)\begin{equation} u = \sum^{\infty}_{j=0} a_j (x,t) \phi_j (z),\quad \theta = 1 + \sum^{\infty}_{j=0} b_j (x,t) \phi_j (z), \end{equation}

where $\{ a_j , b_j \}$ denote unknown expansion coefficients and $\phi _j = ( z - s)^{j+1}$ are basis functions selected to automatically satisfy the substrate boundary conditions (2.14); thus, expansions (2.17) represent power series in $z$ centred around the substrate position, $s ( x )$. An expansion of the vertical flow velocity is unnecessary since an expression for $w$ in terms of $u$ can be found by expressing the continuity equation (2.13) in the form

(2.18)\begin{equation} w ={-} \int^{z}_{s} \frac{ \partial u }{ \partial x }\, \textrm{d}z. \end{equation}

An expansion of the fluid pressure is likewise redundant since an algebraic expression – (2.16) – exists.

Representing $u$ and $\theta$ as power series in $z$ enables a reduction in dimensionality because substituting the expansions (2.17) into the governing equations (2.13) converts each equation into an infinite sum of $\phi _j$; this can be represented mathematically by

(2.19)\begin{equation} R_{({-}1)} (x,t) + R_{0} (x,t) \phi_0 ( \hat{z} ) + R_{1} (x,t) \phi_1 ( \hat{z} ) + \cdots+ R_{J} (x,t) \phi_J ( \hat{z} ) + \cdots= 0 ,\end{equation}

in which the $\phi _j$ depend on the reduced variable $\hat {z} = ( z - s )$ whilst the coefficients $R_j$ are functions of $( x , t )$. A reduction in dimensionality follows because: (i) (2.19) must be satisfied across the whole domain of $\varOmega = \{ x \in [ 0 , 1 ] , \ \hat {z} \in [ 0 , h ] \}$; and (ii) $\{ 1 , \phi _j \}$ constitute the monomial basis of $\hat {z}$ and are therefore linearly independent. Consequently, (2.19) can only be satisfied if $R_j (x,t) = 0$ for $j = -1 , 0 , 1 ,\ldots, J ,\ldots, \infty$. This leads to a reduced equation set made up of the residuals $R_j (x,t)$ for $j = -1 , 0 , 1 ,\ldots, J ,\ldots, \infty$ and the free-surface boundary conditions (2.15) which are functions of $( x , t )$ but independent of $z$. The reduced equations are satisfied by solutions to the expansion coefficients $\{ a_j , b_j \}$.

Needless to say, solving an infinite number of residual equations would be impractical, and so truncations of expansions (2.17) are required; thus, the solutions to $u$ and $\theta$ sought via this approach are approximate rather than exact. Nevertheless, it is at least possible to find asymptotically equivalent solutions within the framework of a long-wave expansion. The smallness of derivatives with respect to $( x , t )$ means approximating such terms using truncations of expansions (2.17) leads to asymptotic solutions which are equivalent to the exact solutions within a bounded parameter space. The truncations needed to approximate space–time derivatives are found from a gradient expansion of the fluid velocity and temperature, taking the form $( u , \theta ) = ( u_0 , \theta _0 ) + {\epsilon } ( u_1 , \theta _1 ) + {\epsilon }^2 ( u_2 , \theta _2 ) + \cdots$.

The first step is to find $( u_0 , \theta _0)$. A potential starting point are the analytical solutions for the case of ‘a flat film flowing down a planar, uniformly heated incline’, found by setting ${\epsilon } = 0$ in (2.13)–(2.15), and attributed to Nusselt (Reference Nusselt1916), who first found such an expression for the velocity profile through an isothermal flat film; they read

(2.20a,b)\begin{equation} u_{Nu} = ( z - s ) ( s + 2 h - z ),\quad \theta_{Nu} = 1 + \left( \frac{1}{1 + Bih} - 1 \right) \frac{ ( z - s ) }{ h }. \end{equation}

Equations (2.20) represent exact solutions which are only strictly valid when ${\epsilon } = 0$; indeed, setting $( u_0 , \theta _0 ) = ( u_{Nu} , \theta _{Nu})$ leads to a Benney expansion in terms of the film thickness $( h )$. The derivation of a reduced asymptotic model which accurately resolves the film dynamics beyond the trivial case requires expansions whose forms are invariant with respect to ${\epsilon }$ – consider the self-similarity function in a Blasius boundary layer (Blasius Reference Blasius1908) – and necessitates a multi-modal parametrisation of the reduced equations (Roberts Reference Roberts1996).

A better starting point is to identify which of the expansion coefficients $\{ a_j ,b_j \}$ can enter at leading order in the gradient expansion, and use these expansion coefficients $\{ a_j , b_j \}$ to represent the amplitudes of the primary flow and heat transfer, i.e. $( u_0 , \theta _0 )$. This is accomplished by deriving a set of conditions from the governing equations for the heated film case which isolate the individual expansion coefficients so their order may be inferred from other terms present in the equations (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000). Since expansions (2.17) are power series in $z$, $\{ a_j , b_j \}$ can be isolated by differentiating equations (2.13b,d) $n - 2$ times with respect to $z$ and only with respect to $z$; this gives

(2.21a)$$\begin{gather} {\epsilon} Re \frac{ \partial^{n-2} }{ \partial z^{n-2} } \left[ u \frac{ \partial u }{ \partial x } + w \frac{ \partial u }{ \partial z } + \frac{ \partial u }{ \partial t }\right] = \frac{ \partial^n u }{ \partial z^n }, \end{gather}$$

(2.21b)$$\begin{gather}{\epsilon} Re\, Pr \frac{ \partial^{n-2} }{ \partial z^{n-2} } \left[ u \frac{ \partial \theta }{ \partial x } + w \frac{ \partial \theta }{ \partial z } + \frac{ \partial \theta }{ \partial t } \right] = \frac{ \partial^n \theta }{ \partial z^n }, \end{gather}$$

where it has been assumed $Re \sim Pr \sim {O} (1)$ and all terms of ${\sim}{O} ( {\epsilon }^2 )$ and higher have been discarded. At this point it is only necessary to determine which $\{ a_j , b_j \}$ are of at least first order in the long-wave expansion. Considering all powers of $z$, the infinite expansions (2.17) yield the following conditions when substituted into (2.21):

(2.22) \begin{align} &\boldsymbol{R} \left[\sum^{n-2}_{k=0} \left( \sum^{\infty}_{j=0} a_j \phi_{j-n+k+2} \prod^{n-k-2}_{i=0} \frac{ ( \,j - i + 2 ) }{ ( \,j + 2 ) } \right) \left( \sum^{\infty}_{j=0} \frac{ \partial }{ \partial x } [ \boldsymbol{r}_j \phi_{j-k} ] \prod^{k}_{i=0} \frac{ ( \,j - i + 2 ) }{ ( \,j + 2 ) } \right) \right. \nonumber\\ &\quad - \sum^{n-2}_{k=0}\left(\sum^{\infty}_{j=0} \frac{ \dfrac{ \partial }{ \partial x } [ a_j \phi_{j-n+k+3} ] }{ ( \,j + 2 ) } \prod^{n-k-2}_{i=0} \frac{ ( \,j - i + 3 ) }{ ( \,j + 3 ) }\right) \left( \sum^{\infty}_{j=0} \boldsymbol{r}_j \phi_{j-k-1} \prod^{k+1}_{i=0} \frac{ ( \,j - i + 2 ) }{ ( \,j + 2 ) } \right) \nonumber\\ &\quad \left.+ \sum^{\infty}_{j=0} \frac{ \partial \boldsymbol{r}_j }{ \partial t } \phi_{j-n+2} \prod^{n-2}_{i=0} \frac{ ( \,j - i + 2 ) }{ ( \,j + 2 ) } \right] = \sum^{\infty}_{j=0} \boldsymbol{r}_j \phi_{j-n} \prod^{n}_{i=0} \frac{ ( \,j - i + 2 ) }{ ( \,j + 2 ) }, \end{align}

for (i) $\boldsymbol {r}_j = a_j$ with $\boldsymbol {R} = {\epsilon } Re$, and (ii) $\boldsymbol {r}_j = b_j$ with $\boldsymbol {R} = {\epsilon } Re\, Pr$.

The right-hand side of (2.22) is proportional to $\{ a_{n-1} , b_{n-1} \}$ because the repeated product only leaves terms of $\{ a_j , b_j \}$ satisfying the inequality $j - n + 2 > 0$, simplifying to $j \geq n - 1$ as $\{ \,j , n \} \in \mathbb {Z}^+$. Meanwhile, the left-hand side of (2.22) is of at least first order in the long-wave expansion due to the presence of space–time derivatives. Stipulating that both sides of (2.22) be of equivalent orders requires $\{ a_{n-1} , b_{n-1} \}$ to be of at least first order. Equations (2.21)–(2.22) hold for $n > 2$: constants of integration appear when $n \leq 2$; implying all $\{ a_j , b_j \}$ for $j \geq 2$ are of at least first-order; leaving only $\{ a_0 , a_1 , b_0 , b_1 \}$ to be candidates of order unity, ${\sim}{O} ( 1 )$, and to represent the amplitudes of the primary flow and heat transfer, i.e. $( u_0 , \theta _0 )$.

On the conjecture that $\{ a_0 , a_1 , b_0 , b_1 \}$ are the leading expansion coefficients, attention is directed to expressing $\{ a_0 , a_1 , b_0 , b_1 \}$ in terms of physically measurable parameters and making sure the power series (2.17) satisfy all of the boundary conditions (2.14)–(2.15). Consequently, a variable transformation of the following form is sought:

(2.23)\begin{equation} (u , w , \theta) \to ( q , h , \vartheta ), \end{equation}

where $q$ is the streamwise flow rate, $h$ is the film thickness and $\vartheta$ is the free-surface temperature. The flow rate $( q )$ is a controllable physical quantity, whilst the film thickness $( h )$ and free-surface temperature $( \vartheta )$ already appear in the free-surface boundary conditions (2.15). The variable transformation can be easily accomplished by finding the expressions which relate $\{ a_j , b_j \}$ to $( q , h , \vartheta )$; in addition to providing physical context, the resulting expressions describe how the velocity expansion coefficients $\{ a_j \}$ and the temperature expansion coefficients $\{ b_j \}$ linearly depend upon one another, respectively.

The fluid velocity and temperature are related to the streamwise flow rate and free-surface temperature, respectively, by

(2.24a,b)\begin{equation} q = \int^{f}_{s} u \, \textrm{d}z, \quad \vartheta = \theta |_{z=f}. \end{equation}

Substituting expansions (2.17) into definitions (2.24) leads to the following relationships:

(2.25a,b)\begin{equation} q = \sum^{\infty}_{j=0} \frac{ h^{j+2} }{ ( \,j + 2) } a_j, \quad \vartheta - 1 = \sum^{\infty}_{j=0} h^{j+1} b_j, \end{equation}

which form the first set of conditions relating $\{ a_j , b_j \}$ to $( q , h , \vartheta )$.

The decision to centre the power series about $s ( x )$ ensures the substrate boundary conditions (2.14) are automatically satisfied by expansions (2.17); however, the free-surface boundary conditions (2.15) are not. Substituting expressions (2.17)–(2.18) into (2.15) yields supplementary conditions which complement those in (2.25). Approximating the kinematic condition (2.15a) using (2.17)–(2.18), and simplifying via condition (2.25), leads to the integral form of the kinematic condition

(2.26)\begin{equation} \frac{ \partial h }{ \partial t } + \frac{ \partial q }{ \partial x } = 0, \end{equation}

describing the relationship between film thickness and flow rate. Equation (2.26) is universal for non-volatile, incompressible film flow and can be derived by integrating the continuity equation (2.13a) between $z = s(x)$ and $z = f(x,t)$, with the lower and upper bounds of the integration given by (2.14a) and (2.15a), respectively, followed by simplifying via the Leibniz integral rule and streamwise flow rate definition.

Since both the $z$-momentum equation (2.13c) and the normal stress condition (2.15b) are eliminated from the equation set (2.13)–(2.15) following a substitution of the fluid pressure – (2.16) – into the $x$-momentum – (2.13); only the shear stress (2.15c) and heat flux (2.15d) boundary conditions remain. With reference to the variable transformation of (2.23), boundary conditions (2.15c,d) can be re-written as

(2.27a)$$\begin{gather} \left. \frac{ \partial u }{ \partial z } \right|_{z=f} = \frac{ 2 {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \left.\dfrac{ \partial u }{ \partial x } \right|_{z=f} }{ 1 - {\epsilon}^2 g } + \dfrac{ {\epsilon}^2 \left[ \dfrac{ \partial^2 q }{ \partial x^2 } - u|_{z=f} \dfrac{ \partial^2 f }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g } - \dfrac{ {\epsilon} \dfrac{Ma}{Ca} \dfrac{ \partial \vartheta }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } [ 1 - {\epsilon}^2 g ] }, \end{gather}$$

(2.27b)$$\begin{gather}\left. \frac{ \partial \theta }{ \partial z } \right|_{z=f} = \frac{ {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } - \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g} }, \end{gather}$$

utilising $\partial w / \partial z = - \partial u / \partial x$ via the continuity equation (2.13a), $w|_{z=f}$ from the kinematic condition (2.15a) and the total spatial derivatives defined along the free surface, namely

(2.28a)$$\begin{gather} \frac{ \partial \vartheta }{ \partial x } = \left. \frac{ \partial \theta }{ \partial x } \right|_{z=f} + \frac{ \partial f }{ \partial x } \left. \frac{ \partial \theta }{ \partial z } \right|_{z=f,} \end{gather}$$

(2.28b)$$\begin{gather}\frac{ \partial }{ \partial x } ( u|_{z=f} ) = \left. \frac{ \partial u }{ \partial x } \right|_{z=f} + \frac{ \partial f }{ \partial x } \left. \frac{ \partial u }{ \partial z } \right|_{z=f,} \end{gather}$$

(2.28c)$$\begin{gather}\frac{ \partial }{ \partial x } ( w|_{z=f} ) = \left. \frac{ \partial w }{ \partial x } \right|_{z=f} + \frac{ \partial f }{ \partial x } \left. \frac{ \partial w }{ \partial z } \right|_{z=f}. \end{gather}$$

Substituting expansions (2.17) into the left-hand side of (2.27a,b) (within the framework of the long-wave expansion, space–time derivatives on the right-hand side of (2.27) are approximated via recursion using truncations of expansions (2.17)) yields

(2.29a)\begin{align} \sum^{\infty}_{j=0} ( \,j + 1 ) h^j a_j &= \frac{ 2 {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \left. \dfrac{ \partial u }{ \partial x } \right|_{z=f} }{ 1 - {\epsilon}^2 g } + \dfrac{ {\epsilon}^2 \left[ \dfrac{ \partial^2 q }{ \partial x^2 } - u|_{z=f} \dfrac{ \partial^2 f }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g }\nonumber\\ &\quad- \frac{ {\epsilon} \dfrac{Ma}{Ca} \dfrac{ \partial \vartheta }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } [ 1 - {\epsilon}^2 g ]}, \end{align}

(2.29b)\begin{align} \sum^{\infty}_{j=0} ( \,j + 1 ) h^j b_j &= \frac{ {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } - \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } }, \end{align}

respectively, which form the second set of conditions relating $\{ a_j , b_j \}$ to $( q , h , \vartheta )$.

Thus far, the presumption from (2.22) is that $\{ a_0 , a_1 , b_0 , b_1 \}$ are the leading expansion coefficients, and four relationships linking $\{ a_j , b_j \}$ to $( q , h ,\vartheta )$ have been derived – (2.25) and (2.29); with an additional condition linking the film thickness and flow rate in (2.26). Therefore, solving for $\{ a_0 , a_1 , b_0 , b_1 \}$ in (2.25) and (2.29) returns expressions for them in terms of $( q , h , \vartheta, a_j , b_j )$ for $j \geq 2$, which read

(2.30a)\begin{align} a_0 &=\frac{ 3 q }{ 2h^3 } (2h) + \frac{1}{2} \sum^{\infty}_{j=2} \frac{ ( \,j + 4 ) ( \,j -1 ) }{ ( \,j + 2 ) } h^j a_j \nonumber\\ &\quad - \frac{1}{2} \left[ \frac{ \left. 2 {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial u }{ \partial x } \right|_{z=f} - \dfrac{ {\epsilon} \frac{Ma}{Ca} \frac{ \partial \vartheta }{\partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } } }{ 1 - {\epsilon}^2 g } + \frac{ {\epsilon}^2 \left[ \dfrac{ \partial^2 q }{ \partial x^2 } - u|_{z=f} \dfrac{ \partial^2 f }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g } \right], \end{align}

(2.30b)\begin{align} a_1 &={-} \frac{3 q}{2h^3} - \frac{3}{4h} \sum^{\infty}_{j=2} \frac{ j ( \,j + 3 ) }{ ( \,j + 2 ) } h^j a_j \nonumber\\ &\quad + \frac{3}{4h} \left[ \frac{ \left. 2 {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial u }{ \partial x } \right|_{z=f} - \dfrac{ {\epsilon} \frac{Ma}{Ca} \frac{ \partial \vartheta }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } } }{ 1 - {\epsilon}^2 g } + \frac{ {\epsilon}^2 \left[ \dfrac{ \partial^2 q }{ \partial x^2 } - u|_{z=f} \dfrac{ \partial^2 f }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g } \right], \end{align}

(2.30c)$$\begin{align} b_0 &= \frac{ ( \vartheta - 1 ) }{ h } + \left[ \frac{ ( \vartheta - 1 ) }{ h } - \frac{ {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } + \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] + \sum^{\infty}_{j=2} ( \,j - 1 ) h^j b_j, \end{align}$$

(2.30d)$$\begin{align}b_1 &={-} \frac{1}{h} \left[ \frac{ \left( \vartheta - 1 \right) }{ h } - \frac{ {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } + \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] - \frac{1}{h} \sum^{\infty}_{j=2} j h^j b_j, \end{align}$$

where $( q , h , \vartheta, Ma, Ca, Bi ) \sim {O} (1)$ and $\{ a_j , b_j \} \sim {O} ( {\epsilon }^n )$ for $j \geq 2$ and $n\geq 1$.

This is advantageous because: (i) replacing $\{ a_0 , a_1 , b_0 , b_1 \}$ in the expansions of $( u ,\theta )$ – (2.17) – with the identities above allows the power series to satisfy not only the substrate boundary conditions but also to algebraically satisfy those at the free surface – leaving only the $x$-momentum (2.13b) and energy (2.13d) equations to be satisfied; and (ii) the leading truncations of the gradient expansion, $( u_0 , \theta _0 )$, can now be found. The identities of $\{ a_0 , a_1 \}$ are asymptotically equivalent, in the long-wave limit, to the amplitude of the parabolic velocity profile derived by Ruyer-Quil et al. (Reference Ruyer-Quil, Scheid, Kalliadasis, Velarde and Zeytounian2005); meanwhile, the identities of $\{ b_0 , b_1 \}$ can be compared with the temperature expansions of Trevelyan et al. (Reference Trevelyan, Scheid, Ruyer-Quil and Kalliadasis2007) and Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020). In the work of Trevelyan et al. (Reference Trevelyan, Scheid, Ruyer-Quil and Kalliadasis2007), the expansion coefficients were not assigned any order with respect to ${\epsilon }$ but an emphasis was placed on the cubic temperature coefficient to ease comparison with a Benney expansion of the fluid temperature; accordingly, a temperature expansion with the same form as (5.10b) in Trevelyan et al. (Reference Trevelyan, Scheid, Ruyer-Quil and Kalliadasis2007) can be retrieved by solving for $b_2$ instead of $b_1$ in (2.30c,d). Unfortunately, the single-variable heat transfer model proposed by Trevelyan et al. (Reference Trevelyan, Scheid, Ruyer-Quil and Kalliadasis2007) is inconsistent with the long-wave expansion because it was derived following the Galerkin method; a consistent evolution equation for $\vartheta$ can be obtained by following either the current approach or that of Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020). In a similar fashion, a temperature expansion with the same form as Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020) can be retrieved by introducing the evaluation of the energy equation (2.13c) at the free surface as a supplementary boundary condition; the latter can be written as

(2.31)\begin{align} \left. \frac{ \partial^2 \theta }{ \partial z^2 } \right|_{z=f} = \frac{ {\epsilon} Re\,Pr \left[ \dfrac{ \partial \vartheta }{ \partial t } + u|_{z=f} \dfrac{ \partial \vartheta }{ \partial x } \right] + {\epsilon}^2 \left[ \dfrac{ \partial^2 f }{ \partial x^2 } \left. \dfrac{ \partial \theta }{ \partial z } \right|_{z=f} + 2 \dfrac{ \partial f }{ \partial x } \dfrac{ \textrm{d} }{ {\textrm{d}\kern0.7pt x} } \left( \left. \dfrac{ \partial \theta }{ \partial z } \right|_{z=f} \right) - \dfrac{ \partial^2 \vartheta }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g }. \end{align}

Substituting the power series (2.17b) into the left-hand side of (2.31), with $\{ b_0 , b_1 \}$ given by (2.30c,d), and solving for $b_2$ returns a temperature expansion with the same form as Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020). Equation (2.31) is left out of the current analysis because the space–time derivatives composing the right-hand side indicate it is purely a first-order contribution; the present objective is to find $( u_0 , \theta _0 )$.

To write the expansions at leading order, expressions (2.30) are substituted into the power series (2.17) and all terms of order ${\sim}{O}( {\epsilon })$ and higher are discarded; this yields

(2.32a)$$\begin{gather} u_0 = \frac{3 q}{2h^3} ( z - s ) ( s + 2h - z ), \end{gather}$$

(2.32b)$$\begin{gather}\theta_{para} = 1 + \frac{ ( \vartheta - 1 ) }{ h } ( z - s ) + \frac{1}{h} \left[ \frac{ ( \vartheta - 1 ) }{ h } + \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] ( z - s ) ( s + h - z ), \end{gather}$$

in which surface curvature $( {\epsilon }^2 g )$ has been retained, so its effect on the heat transfer through the free surface is accurately captured, but thermo-capillarity has been excluded as the Marangoni effect is purely a deviant mechanism in uniformly heated film flow. The leading temperature expansion (2.32b) is labelled $\theta _{para}$ because it is parabolic and thus a departure from the linear temperature expansion used by Kalliadasis et al. (Reference Kalliadasis, Demekhin, Ruyer-Quil and Velarde2003a,Reference Kalliadasis, Kiyashko and Demekhinb), Ruyer-Quil et al. (Reference Ruyer-Quil, Scheid, Kalliadasis, Velarde and Zeytounian2005), D'Alessio et al. (Reference D'Alessio, Pascal, Jasmine and Ogden2010) and Sterman-Cohen & Oron (Reference Sterman-Cohen and Oron2020).

A long-standing assumption in the modelling of heated film flow has been that the leading approximation to the temperature expansion must be linear because the Nusselt temperature distribution – (2.20b) – is linear. This line of thinking would demand the parabolic correction in (2.32b) be of first order with respect to ${\epsilon }$ and reduce the leading temperature expansion to the linear form offered by Ruyer-Quil et al. (Reference Ruyer-Quil, Scheid, Kalliadasis, Velarde and Zeytounian2005). There is merit in this assumption because the parabolic coefficient $b_1$ does tend to zero in the case of a ‘flat film flowing down a planar, uniformly heated incline’ (due to the free-surface temperature approaching its flat-film solution), or expressed mathematically

(2.33a,b)\begin{equation} \lim_{{\epsilon} \to 0 } \vartheta = \frac{ 1 }{ 1 + Bi h } + {O} ( {\epsilon} ), \quad \lim_{{\epsilon} \to 0} \left[ \frac{ ( \vartheta - 1 ) }{ h } + \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] + {O} ( {\epsilon} ) = 0. \end{equation}

Nevertheless, there are two issues with this assumption: (i) the parabolic correction in (2.32b) contains only quantities which are of order unity, $( h , \vartheta, Bi ) \sim {O} (1)$; and (ii) a linear temperature approximation at leading order produces inaccurate predictions except when the free-surface temperature $( \vartheta )$ remains close to its flat-film solution. To resolve this quandary, it must be understood why the modelling approach of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000) is superior to the traditional approach of Benney (Reference Benney1966).

The superiority of the current modelling approach stems from the fact that expressions (2.32) are approximations to the general solutions of the fluid velocity and temperature at leading order; in contrast, the Nusselt solutions – (2.20) – are particular solutions corresponding to the special case of ${\epsilon } = 0$. Consequently, expressions (2.32) are applicable at any value of ${\epsilon }$ whereas (2.20) are strictly only valid when ${\epsilon } = 0$. An important test is to check whether the general solutions – (2.32) – reduce to the particular solutions – (2.20) – when the special case of ${\epsilon } = 0$ is met.

Writing the momentum (2.13b) and energy (2.13c) equations with ${\epsilon } = 0$, returns

(2.34a,b)\begin{equation} \frac{ \partial^2 u }{ \partial z^2 } ={-}2,\quad \frac{ \partial^2 \theta }{ \partial z^2 } = 0, \end{equation}

which when evaluated using expressions (2.32), with ${\epsilon }^2 g = 0$; finds $q$ and $\vartheta$ to equal

(2.35a,b)\begin{equation} q_{Nu} = \frac{2}{3} h^3, \quad \vartheta_{Nu} = \frac{ 1 }{ 1 + Bi h }. \end{equation}

Expressions (2.35) are the flat-film (or Nusselt) solutions to the flow rate and free-surface temperature, respectively; substituting these solutions into expressions (2.32), with ${\epsilon }^2 g = 0$, reduces the leading expansions to the Nusselt solutions – (2.20) – as required.

Of course, the trivial case $( {\epsilon } = 0 )$ is not the point of interest; the aim is to find particular solutions for when ${\epsilon } \neq 0$. The general solutions are applicable at any value of ${\epsilon }$ and can therefore be used to transform the $x$-momentum (2.13b) and energy (2.13d) equations into an infinite set of residual equations which are of reduced dimensionality and which automatically satisfy the boundary conditions (2.14)–(2.15). The infinite set of residuals can then be truncated via a gradient expansion, $( u , \theta ) = ( u_0 , \theta _0 ) + {\epsilon } ( u_1 , \theta _1 ) + {\epsilon }^2 ( u_2 , \theta _2 ) + \cdots$ where ${\epsilon } < 1$, to a finite set whose solutions are asymptotically equivalent to those of the full equation set (2.13)–(2.15). In contrast, the Nusselt solutions are only valid when ${\epsilon } = 0$ and proceeding to set $( u_0 , \theta _0 ) = ( u_{Nu} , \theta _{Nu} )$ is tantamount to assuming the powers of ${\epsilon }$ are linearly independent in the $x$-momentum (2.13b) and energy (2.13d) equations since it leads to each of the coefficients of ${\epsilon }^n$ being equal to zero, i.e. $\{ 0 \} + {\epsilon } \{ 0 \} + \cdots = 0$. Assuming the powers of ${\epsilon }$ to be linearly independent is in direct contradiction to the gradient expansion which approximates the fluid velocity and temperature through a summation of the powers of ${\epsilon }$ and thus relies upon them being linearly dependent. In fact, if the powers of ${\epsilon }$ are treated as being linearly independent in the gradient expansion then this would necessitate $( u_n , \theta _n ) \equiv 0$ for $n \geq 1$ and $( u , \theta ) = ( u_{Nu} , \theta _{Nu} )$; equivalent to setting ${\epsilon } = 0$ and providing some explanation as to why the Benney expansion is only consistent in a narrow neighbourhood about the Nusselt solutions.

Obviously, the powers of ${\epsilon }$ must always be treated as linearly dependent; the leading expansions – expressions (2.32) – allow for this whereas the Nusselt solutions – (2.20) – do not. Interestingly, this detail facilitates the existence of leading-order contributions which are balanced entirely against a summation of higher-order terms; this is completely in keeping with the long-wave expansion which only assumes space–time derivatives are small, thus a large contribution to the film dynamics can be equal to the sum of many smaller contributions; such a leading-order contribution would vanish in the limit of ${\epsilon } \to 0$. The concept of evanescent leading-order contributions has been touched upon previously in the literature: Ooshida (Reference Ooshida1999) theorised the existence of a drag–inertia regime which makes a large contribution to the velocity field but only when the fluid inertia is non-negligible – to explain why the Nusselt solution is an inadequate description of the flow dynamics beyond the trivial case; and Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020) introduced nonlinear relaxation modes to the temperature expansion designed to satisfy Newton's law of cooling in the long-wave limit whilst simultaneously vanishing in this same limit. In either case, retaining these evanescent leading-order contributions has been shown to drastically increase the predictive capability of the models in question.

The discussion above explains how the parabolic correction in expression (2.32b) can be of leading order with respect to ${\epsilon }$ despite vanishing in the limit of ${\epsilon } \to 0$; to showcase this behaviour the temperature predictions of two models are compared in § 3, the first assumes $b_1 \sim {O} ( {\epsilon } )$ whilst the second allows $b_1 \sim {O} (1)$. Perhaps the most paradoxical behaviour of $b_1$ is seen in the problem of film flow down planar inclines, $s(x) = 0$ $\forall x$; where the evaluation of the energy equation at the substrate, $\partial ^2 \theta / \partial z^2 |_{z=s} = 0$, reveals $b_1 = 0$ $\forall {\epsilon }$. However, it is easily seen from (2.30d) that $b_1$ is itself an infinite expansion, made up of terms belonging to different orders of the long-wave expansion. Therefore, at $n$th order, higher-order terms in (2.13c) must be approximated by $\tilde {b}_1$ where $b_1 = \tilde {b}_1 + {O} ( {\epsilon }^n )$, not $b_1$ itself. This means that even when $b_1 = 0$, it is possible for $\tilde {b}_1 \neq 0$. As a matter of fact one should anticipate $\tilde {b}_1 \neq 0$ when ${\epsilon } \neq 0$. Consider (2.33), the square-bracketed term in (2.33b) approaches zero in the limit of ${\epsilon } \to 0$ because $\vartheta$ approaches its flat-film solution in (2.33a). At first-order, $\tilde {b}_1$ is given by the square-bracketed term of (2.33b), thus the only way $\tilde {b}_1$ can equal zero is if $\vartheta = \vartheta _{Nu}$. When ${\epsilon } \neq 0$, $\vartheta \neq \vartheta _{Nu}$ which inevitably means $\tilde {b}_1 \neq 0$.

The leading expansions – expressions (2.32) – were derived on that basis that $\{ a_0 , a_1 , b_0 , b_1 \}$ are the leading expansion coefficients, however, since the powers of ${\epsilon }$ are linearly dependent, it is possible formulate leading-order approximations from any permutation of the expansion coefficients $\{ a_j , b_j \}$. Accordingly, expression (2.32b) is the lowest-degree polynomial approximation of the temperature field to leading order and thus reveals the heat transfer through a heated film to be characteristically nonlinear. It is impossible for a linear temperature expansion to encapsulate all of the thermodynamics in heated film flow when ${\epsilon } \neq 0$ because the general solution must satisfy a total of three boundary conditions and a linear approximation can only ever satisfy two. Alternative leading temperature approximations would be the single-variable temperature expansion of Trevelyan et al. (Reference Trevelyan, Scheid, Ruyer-Quil and Kalliadasis2007) – (5.10b) of their paper with $i \leq 1$ – which uses a cubic correction; or the first temperature ansatz proposed by Cellier & Ruyer-Quil (Reference Cellier and Ruyer-Quil2020) – (17) of their paper – which uses a polynomial correction.

With the necessity of the parabolic temperature profile – expression (2.32b) – established, focus is now directed at reducing the dimensionality of the complete governing equations which contain terms up to fourth order. A substitution of the fluid pressure (2.16) and vertical velocity (2.18) expressions, eliminate both the continuity (2.13a) and $z$-momentum (2.13c) equations; leaving only the $x$-momentum (2.13b) and energy (2.13d) equations. Therefore, it is the reduced forms of (2.13b,d) that accompany (2.26) in forming the reduced asymptotic model.

In the ideal scenario, every term involving $( u , \theta )$ would be approximated to ${\sim}{O} ( {\epsilon }^n )$ so the reduced model is asymptotically equivalent to the equation set (2.13)–(2.15) at $n$th order in the long-wave expansion; Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000) refer to this as the complete model. However, a complete model becomes cumbersome at second order as the gradient expansion stipulates that truncations of the expansions (2.17) at ${\sim}{O} ( {\epsilon }^n )$ include all $\{ a_j , b_j \}$ up to $j \leq 4n + 1$. This results in the space–time derivatives of $\{ a_j , b_j \}$ for $j \geq 2$ appearing at ${\sim}{O} ( {\epsilon }^2 )$; the behaviour of these derivatives is unknown making their algebraic elimination impossible. As a result, a complete model at second order, or higher, consists of more than three equations in terms of the physically ambiguous $\{ a_j , b _j \}$ for $j \geq 2$. Naturally, it is desirable to eliminate all $\{ a_j , b_j \}$ via a simplification. The easiest and most straightforward step is to simply discard all the derivatives of $\{ a_j , b_j \}$ for $j \geq 2$, along with the products of $\{ a_j , b_j \}$ for $j \geq 2$ with any other derivatives. This leads to what has become known as the simplified model (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000) and is justifiable, even in the context of undulating substrate and moderate Reynolds number, because the flow remains primarily parabolic and the film dynamics is chiefly described by the leading velocity expansion (Wang Reference Wang1984; Ruyer-Quil et al. Reference Ruyer-Quil, Scheid, Kalliadasis, Velarde and Zeytounian2005).

The fourth-order simplified model presented here is derived by approximating higher-order terms in (2.13b,d) using the leading expansions (2.32). With the leading terms in (2.13b,d) – the double derivatives with respect to $z$ which introduce only linear terms in $\{ a_j , b_j \}$ – obtained using the complete expansions (2.17), this leads to

(2.36a)\begin{align} &{[} A_0 - 2 a_1 ] + [ A_1 - 6 a_2 ] \phi_0 + [ A_2 - 12 a_3 ] \phi_1 + [ A_3 - 20 a_4 ] \phi_2 + [ A_4 - 30 a_5 ] \phi_3 \nonumber\\ &\quad + [ A_5 - 42 a_6 ] \phi_4 + [ A_6 - 56 a_7 ] \phi_5 - \sum^{\infty}_{m=8} [ m ( m + 1 ) a_{m} ]\phi_{m-2} = 0, \end{align}

(2.36b)\begin{align} &{[} B_0 - 2 b_1 ] + [ B_1 - 6 b_2 ] \phi_0 + [ B_2 - 12 b_3 ] \phi_1 + [ B_3 - 20 b_4 ] \phi_2 + [ B_4 - 30 b_5 ] \phi_3\notag\\ &\quad - \sum^{\infty}_{n=6} [ n ( n + 1 ) b_{n} ] \phi_{n-2} = 0; \end{align}

the expressions for $\{ A_m , B_n \}$ are provided in the online supplementary material in terms of $( s , \tilde {a}_0 , \tilde {a}_1 , \tilde {b}_0 , \tilde {b}_1 )$ – with $a_j = \tilde {a}_j + {O} ( {\epsilon } )$ and $b_j = \tilde {b}_j + {O} ( {\epsilon } )$. Setting the coefficients of $\phi _j$ to zero – as per (2.19) – and solving for $\{ a_j , b_j \}$ for $j \geq 2$, returns

(2.37a[i,ii])$$\begin{gather} a_7 = A_6/56,\qquad a_m = 0, \quad \text{for} \ m \geq 8, \end{gather}$$

(2.37b[i,ii])$$\begin{gather}a_6 = A_5/42,\qquad b_n = 0, \quad \text{for}\ n \geq 6, \end{gather}$$

(2.37c[i,ii])$$\begin{gather}a_5 = A_4/30,\qquad b_5 = B_4/30, \end{gather}$$

(2.37d[i,ii])$$\begin{gather}a_4 = A_3/20,\qquad b_4 = B_3/20, \end{gather}$$

(2.37e[i,ii])$$\begin{gather}a_3 = A_2/12,\qquad b_3 = B_2/12, \end{gather}$$

(2.37f[i,ii])$$\begin{gather}a_2 = A_1/6,\qquad b_2 = B_1/6, \end{gather}$$

from which the identities of $\{ a_j , b_j \}$ for $j \geq 2$, valid to first order in the long-wave expansion, can be ascertained in terms of $( s , \tilde {a}_0 , \tilde {a}_1 , \tilde {b}_0 , \tilde {b}_1 )$ – see Appendices A and B. With the coefficients of $\phi _j$ set to zero, (2.36) can be re-written as

(2.38a)$$\begin{gather} \dfrac{ 2h }{3} A_0 + \dfrac{2 q}{h^2} + \sum^{\infty}_{j=2} \dfrac{ j ( \,j + 3 ) }{ ( \,j + 2 ) } h^j a_j = \dfrac{ 3 {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial }{ \partial x } \left[ \dfrac{q}{h} \right] - \dfrac{ {\epsilon} \frac{Ma}{Ca} \frac{ \partial \vartheta }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } } }{ 1 - {\epsilon}^2 g } + \frac{ {\epsilon}^2 \left[ \dfrac{ \partial^2 q }{ \partial x^2 } - \dfrac{3}{2} \dfrac{q}{h} \dfrac{ \partial^2 f }{ \partial x^2 } \right] }{ 1 + {\epsilon}^2 g }, \end{gather}$$

(2.38b)$$\begin{gather}\dfrac{h }{2} B_0 + \frac{ ( \vartheta - 1 ) }{ h } - \frac{ {\epsilon}^2 \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } + \frac{ Bi \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } + \sum^{\infty}_{j=2} j h^j b_j = 0, \end{gather}$$

in which $\{ a_1 , b_1 \}$ have been replaced by their complete identities from (2.30). Substituting $\{ a_j , b_ j \}$ for $j \geq 2$ from Appendices A and B; $\{ A_0 , B_0 \}$ from the supplementary material; and replacing $\{ \tilde {a}_0 , \tilde {a}_1 , \tilde {b}_0 , \tilde {b}_1 \}$ with asymptotically equivalent expressions in terms of $( q , h , \vartheta )$ – found by truncating equations (2.30) at ${\sim}{O} (1)$ – yields

(2.39a)\begin{align} &{\epsilon} Re \left[ \frac{ \partial q }{ \partial t } - \frac{9}{7} \frac{q^2}{h^2} \frac{ \partial h }{ \partial x } + \frac{17}{7} \frac{q}{h} \frac{ \partial q }{ \partial x } \right] + {\epsilon}^3 Re \left[ \frac{29}{56} \frac{ \partial q }{ \partial t } \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{107h}{336} \frac{ \partial q }{ \partial t } \frac{ \partial^2h }{ \partial x^2 } + \frac{73 q^2}{128} \frac{ \textrm{d}^3 s }{ {\textrm{d}\kern0.7pt x}^3 }\right. \nonumber\\ &\quad + \frac{29}{32} \left( \frac{ \partial q }{ \partial x } \right)^2 \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{509}{672} \left( \frac{ \partial q }{ \partial x } \right)^2 \frac{ \partial h }{ \partial x } - \frac{187 q}{1344} \frac{ \partial^2 q }{ \partial x^2 } \frac{ \partial h }{ \partial x } + \frac{277 q}{128} \frac{ \partial q }{ \partial x } \frac{ \textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } - \frac{15 q}{64} \frac{ \partial^2 q }{ \partial x^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \nonumber\\ &\quad + \frac{163 q^2}{448h} \frac{ \partial^2h }{ \partial x^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{185 q^2}{896h} \frac{ \partial h }{ \partial x } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } + \frac{9 q^2}{7h} \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{11h^2}{56} \frac{ \partial^3 q }{ \partial x^2 \partial t } + \frac{17 q}{7h} \frac{ \partial q }{ \partial x } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2\nonumber\\ &\quad + \frac{1205 q^2}{2688} \frac{ \partial^3h }{ \partial x^3 } - \frac{601h q}{1008} \frac{ \partial^3 q }{ \partial x^3 } - \frac{683h}{1008} \frac{ \partial^2 q }{ \partial x^2 } \frac{ \partial q }{ \partial x } + \frac{297 q}{224h} \frac{ \partial q }{ \partial x } \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{4633 q}{2688} \frac{ \partial q }{ \partial x } \frac{ \partial^2h }{ \partial x^2 }\nonumber\\ &\quad - \frac{555 q^2}{896h} \frac{ \partial^2h }{ \partial x^2 } \frac{ \partial h }{ \partial x } + \frac{11}{8} \frac{ \partial q }{ \partial t } \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{ \partial q }{ \partial t } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 + \frac{5h}{12} \frac{ \partial q}{ \partial t } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } + \frac{115 q^2}{672h^2} \left( \frac{ \partial h }{ \partial x } \right)^3 \nonumber\\ &\quad \left.-\, \frac{11h}{56} \frac{ \partial^2 q }{ \partial x \partial t } \frac{ \partial h }{ \partial x } - \frac{163 q^2}{224h^2} \left( \frac{ \partial h }{ \partial x } \right)^2 \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{253 q}{672h} \frac{ \partial q }{ \partial x } \left( \frac{ \partial h }{ \partial x } \right)^2 - \frac{9 q^2}{7h^2} \frac{ \partial h }{ \partial x } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 \right]\nonumber\\ &\quad + {\epsilon}^4 \left[ \frac{15}{4h} \frac{ \partial q }{ \partial x } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^3 + \frac{5 q}{2h^2} \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^4 - \frac{30 q}{7h^2} \left( \frac{ \partial h }{ \partial x } \right)^4 + \frac{61}{14h} \frac{ \partial q }{ \partial x } \left( \frac{ \partial h }{ \partial x } \right)^3 - \frac{q}{4} \frac{ \partial h }{ \partial x } \frac{ \textrm{d}^3 s }{ {\textrm{d}\kern0.7pt x}^3 }\right. \nonumber\\ &\quad - 7 \frac{ \partial q }{ \partial x } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{q}{4} \frac{ \partial^3h }{ \partial x^3 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{5h}{3} \frac{ \partial^2 q }{ \partial x^2 } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } + \frac{49 q}{4h} \frac{ \partial h }{ \partial x } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{107 q h}{336} \frac{ \partial^4h }{ \partial x^4 } \nonumber\\ &\quad - \frac{5 q h}{12} \frac{ \textrm{d}^4 s }{ {\textrm{d}\kern0.7pt x}^4 } - \frac{5h}{3} \frac{ \partial q }{ \partial x } \frac{ \textrm{d}^3 s }{ {\textrm{d}\kern0.7pt x}^3 } + \frac{113 q}{28h} \frac{ \partial^2h }{ \partial x^2 } \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{19 q}{4h} \left( \frac{ \partial h }{ \partial x } \right)^2 \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } - \frac{23}{4} \frac{ \partial^2 q }{ \partial x^2 } \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \nonumber\\ &\quad + \frac{3 q}{7} \frac{ \partial^3h }{ \partial x^3 } \frac{ \partial h }{ \partial x } + \frac{45 q}{8h} \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 + \frac{7 q}{8} \frac{ \partial^2h }{ \partial x^2 } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } - \frac{107h}{84} \frac{ \partial q }{ \partial x } \frac{ \partial^3h }{ \partial x^3 } - \frac{19}{4} \frac{ \partial^2 q }{ \partial x^2 } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 \nonumber\\ &\quad - \frac{13}{7} \frac{ \partial^2 q }{ \partial x^2 } \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{12 q}{h} \frac{ \partial^2h }{ \partial x^2 } \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{29}{2h} \frac{ \partial q }{ \partial x } \left( \frac{ \partial h }{ \partial x } \right)^2 \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{53 q}{4h^2} \left( \frac{ \partial h }{ \partial x } \right)^2 \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 \nonumber\\ &\quad + \frac{53 q}{56} \left( \frac{ \partial^2h }{ \partial x^2 } \right)^2 - \frac{q}{2} \left( \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } \right)^2 - \frac{17}{14} \frac{ \partial q }{ \partial x } \frac{ \partial^2h }{ \partial x^2 } \frac{ \partial h }{ \partial x } + \frac{11h}{28} \frac{ \partial^3 q }{ \partial x^3 } \frac{ \partial h }{ \partial x } - \frac{14 q}{h^2} \left( \frac{ \partial h }{ \partial x } \right)^3 \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \nonumber\\ &\quad - \frac{13}{4} \frac{ \partial q }{ \partial x } \frac{ \partial h }{ \partial x } \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } - \frac{3 q}{2} \frac{ \textrm{d}^3 s }{ {\textrm{d}\kern0.7pt x}^3 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{181h}{168} \frac{ \partial^2 q }{ \partial x^2 } \frac{ \partial^2h }{ \partial x^2 } + \frac{59}{4h} \frac{ \partial q }{ \partial x } \frac{ \partial h }{ \partial x } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 + \frac{11h^2}{56} \frac{ \partial^4 q }{ \partial x^4 } \nonumber\\ &\quad \left.- \frac{13}{4} \frac{ \partial q }{ \partial x }\frac{ \partial^2h }{ \partial x^2 } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{59 q}{8h} \frac{ \partial^2h }{ \partial x^2 } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 - \frac{5 q}{4h^2} \frac{ \partial h }{ \partial x } \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^3 \right] + \frac{ {\epsilon} \dfrac{5}{4} \dfrac{Ma}{Ca} \dfrac{ \partial \vartheta }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } [ 1 - {\epsilon}^2 g ] } \nonumber\\ &\quad + {\epsilon}^2 \left[ \frac{5 q}{h^2} \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 + \frac{15 q}{8h} \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } + \frac{15}{4h} \frac{ \partial q }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{33}{4h} \frac{ \partial q }{ \partial x } \frac{ \partial h }{ \partial x } - \frac{31 q}{4h^2} \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{33 q}{8h} \frac{ \partial^2h }{ \partial x^2 }\right.\nonumber\\ &\quad \left.- \frac{13}{4} \frac{ \partial^2 q }{ \partial x^2 } - \frac{5 q}{4h^2} \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{15}{4} \frac{ \dfrac{ \partial f }{ \partial x } \dfrac{ \partial }{ \partial x } \left[ \dfrac{q}{h} \right] }{ 1 - {\epsilon}^2 g } - \frac{15}{8} \frac{ \dfrac{2}{3} \dfrac{ \partial^2 q }{ \partial x^2 } - \dfrac{q}{h} \dfrac{ \partial^2 f }{ \partial x^2 } }{ 1 + {\epsilon}^2 g } \right] + {\epsilon} \frac{5h}{3} \frac{ \partial f }{ \partial x } \,\cot {\beta} \nonumber\\ &\quad + \frac{5 q}{2h^2} - \frac{5h}{3} - {\epsilon}^3 \frac{5h}{6} \frac{ \partial^2 }{ \partial x^2 } \left[ \frac{ \left( 1 - Ma \vartheta \right) }{ Ca } \frac{ \dfrac{ \partial f }{ \partial x } }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] = 0, \end{align}

(2.39b)\begin{align} &{\epsilon} Re\, Pr \left[ \frac{ Bi h }{ \sqrt{ 1 + {\epsilon}^2 g } } \left( \frac{1}{5} \frac{ \partial \vartheta }{ \partial t } - \frac{11}{50} \frac{ \vartheta }{ h } \frac{ \partial q }{ \partial x } + \frac{6}{25} \frac{ q \vartheta }{ h^2 } \frac{ \partial h }{ \partial x } + \frac{6}{25} \frac{ q }{ h } \frac{ \partial \vartheta }{ \partial x } \right) + \frac{3}{25} \frac{ ( \vartheta - 1 ) }{h} \frac{ \partial q }{ \partial x }\right. \nonumber\\ &\quad \left.+ \frac{33}{25} \frac{q}{h} \frac{ \partial \vartheta }{ \partial x } + \frac{ \partial \vartheta }{ \partial t } \right] + \frac{ {\epsilon}^3 Re\, Pr\, Bi h \vartheta }{ [ 1 + {\epsilon}^2 g ]^{3/2} } \frac{ \partial f }{ \partial x } \left[ \frac{1}{5} \frac{ \partial^2 q }{ \partial x^2 } - \frac{6}{5} \frac{q}{h} \frac{ \partial^2 f }{ \partial x^2 } \right] \nonumber\\ &\quad + {\epsilon}^2 \left[ \frac{4 ( \vartheta - 1 )}{5h} \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 } + \frac{12 ( \vartheta - 1 )}{5h^2} \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 + \frac{2 ( \vartheta - 1 )}{5h^2} \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{2 ( \vartheta - 1 )}{5h} \frac{ \partial^2h }{\partial x^2 }\right. \nonumber\\ &\quad \left.+\, \frac{8 ( \vartheta - 1 )}{5h^2} \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{ \partial^2 \vartheta }{ \partial x^2 } + \frac{8}{5h} \frac{ \partial \vartheta }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{4}{5h} \frac{ \partial h }{ \partial x } \frac{ \partial \vartheta }{ \partial x } - \frac{12}{5h} \frac{ \dfrac{ \partial f }{ \partial x } \dfrac{ \partial \vartheta }{ \partial x } }{ 1 + {\epsilon}^2 g } \right] \nonumber\\ &\quad + {\epsilon}^2 \frac{ Bi }{ \sqrt{ 1 + {\epsilon}^2 g } } \left[ \frac{6 \vartheta}{5h} \left( \frac{ \partial h }{ \partial x } \right)^2 + \frac{16 \vartheta}{5h} \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } + \frac{12 \vartheta}{5h} \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 - \frac{3 \vartheta}{5} \frac{ \partial^2h }{ \partial x^2 } - \frac{2 \vartheta}{5} \frac{\textrm{d}^2 s }{ {\textrm{d}\kern0.7pt x}^2 }\right. \nonumber\\ &\quad \left.- \frac{6}{5} \frac{ \partial h }{ \partial x } \frac{ \partial \vartheta }{ \partial x } - \frac{4}{5} \frac{ \partial \vartheta }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } - \frac{h}{5} \frac{ \partial^2 \vartheta }{ \partial x^2 } \right] + {\epsilon}^4 \frac{ Bi }{ [ 1 + {\epsilon}^2 g ]^{3/2} } \left[ \frac{h \vartheta}{5} \frac{ \partial f }{ \partial x } \frac{ \partial^3 f }{ \partial x^3 }\right. \nonumber\\ &\quad \left.+ \frac{6 \vartheta}{5} \left( \frac{ \partial h }{ \partial x } \right)^2 \frac{ \partial^2 f }{\partial x^2 } + 2 \vartheta \frac{ \partial h }{ \partial x } \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} }\frac{ \partial^2 f }{ \partial x^2 } + \frac{4 \vartheta}{5} \left( \frac{\textrm{d} s }{ {\textrm{d}\kern0.7pt x} } \right)^2 \frac{ \partial^2 f }{ \partial x^2 } +\frac{2h}{5} \frac{ \partial f }{ \partial x } \frac{ \partial^2 f }{ \partial x^2 } \frac{ \partial \vartheta }{ \partial x } \right] \nonumber\\ &\quad + {\epsilon}^4 \frac{2 Bi h \vartheta}{5} \frac{ 1 - {\epsilon}^2 g }{ [ 1 + {\epsilon}^2 g ]^{5/2} } \left( \frac{ \partial^2 f }{ \partial x^2 } \right)^2 + \frac{12}{5h^2} \left[ \vartheta - 1 + \frac{ Bi h \vartheta }{ \sqrt{ 1 + {\epsilon}^2 g } } \right] = 0, \end{align}

which together with (2.26) constitutes the fourth-order simplified model utilised in § 3. To be clear, the prefix simplified used here refers to the fact that space–time derivatives of $\{ a_j , b_j \}$ for $j \geq 2$ have been neglected in the above formulation.

Accordingly, the simplified model is only asymptotically equivalent with the full equation set – (2.13)–(2.15) – to first order in the long-wave expansion. From a physical standpoint, the expansion coefficients $\{ a_j , b_j \}$ for $j \geq 2$ describe deviations of the fluid velocity and temperature away from the primary flow and heat transfer. As such, these deviations are retained in the simplified formulation; however, their spatial and temporal evolution is lost. This makes the simplified formulation well suited for analysing film stability down flat plates when any deviation from the primary parabolic flow signals instability; however, beyond first order in the long-wave expansion, inaccuracies can be expected. Even so, it is known the flow rate, film thickness and free-surface temperature still dictate the flow dynamics at moderate Reynolds number and so such inaccuracies are often small and inconsequential (Ruyer-Quil et al. Reference Ruyer-Quil, Scheid, Kalliadasis, Velarde and Zeytounian2005). In the particular case of film flow down wavy substrate, retaining terms up to fourth order makes the present formulation applicable to problems involving large substrate amplitude when the effects of vertical inertia and viscosity are significant – provided the ratio of the substrate amplitude to the substrate wavelength remains small enough for the primary flow to remain parabolic, i.e. $A/L \sim {O} ( {\epsilon } )$. As in the isothermal flow case (Veremieiev & Wacks Reference Veremieiev and Wacks2019), but without formally establishing the necessity of a parabolic temperature profile, (2.39) can be arrived at using weighted residuals. The procedure for this being: (i) higher-order terms involving $( u , \theta )$ are approximated by the leading expansions (2.32); (ii) the double derivatives with respect to $z$ are evaluated using the complete expansions (2.17) with $\{ a_0 , a_1 , b_0 , b_1 \}$ given by (2.30); (iii) weighting of the momentum equation (2.13b) by the Nusselt parabolic velocity profile, $\tilde {w}_{u} = ( z - s ) ( s + 2h - z )$, and of the energy equation (2.13c) by the Nusselt linear temperature profile, $\tilde {w}_{\theta } = ( z - s )$; and finally (iv) integrating the equations between $z = s(x)$ and $z = s(x) + h(x,t)$.

To distinguish the underpinning of the current analysis, (2.26) and (2.39) are referred to henceforth as the reduced asymptotic model (RAM) based on a parabolic temperature profile through the film, or RAM–$\theta _{para}$ for short.

2.3. Linear stability analysis

An essential requirement of the proposed RAM–$\theta _{para}$ is that it returns the critical Reynolds number for the case of a ‘flat film flowing down a planar, uniformly heated incline.’ The accompanying analysis, see Appendix C, leads to the following expression:

(2.40)\begin{equation} Re - \frac{5}{4} \cot \beta + \frac{15}{16} \frac{Ma}{Ca} \frac{ Bi }{ ( 1 + Bi )^2 } = 0, \end{equation}

equivalent to that obtained by Goussis & Kelly (Reference Goussis and Kelly1991) and who were the first to do so. Expression (2.40) is derived in the long-wave limit; thus, it is only valid when the wavenumber of the most unstable mode tends to zero as $Re \to 0$ (Kalliadasis et al. Reference Kalliadasis, Demekhin, Ruyer-Quil and Velarde2003a). It is similarly recovered by assuming the leading temperature expansion is linear – see Scheid et al. (Reference Scheid, Ruyer-Quil, Kalliadasis, Velarde and Zeytounian2005) – because nonlinear behaviour is neglected in the linearised theory.

The stability of film flow down heated wavy substrate is explored using the Floquet theory approach of Trifonov (Reference Trifonov2014), extended to the thermal problem of interest here. Perturbing and linearising the RAM–$\theta _{para}$ (2.26) and (2.39), yields

(2.41a)$$\begin{gather} \frac{ \partial \hat{h} }{ \partial t } + \frac{ \partial \hat{q} }{ \partial x } = 0, \end{gather}$$

(2.41b)$$\begin{gather}\sum^{k=2}_{k=0} \alpha_k (x) \frac{ \partial^{k+1} \hat{q} }{ \partial t \partial x^{k} } + \sum^{k=4}_{k=0} \left[ \beta_k (x) \frac{ \partial^{k} \hat{h} }{ \partial x^{k} } + \gamma_k (x) \frac{ \partial^{k} \hat{q} }{ \partial x^{k} } + \xi_k (x) \frac{ \partial^{k} \hat{\vartheta} }{ \partial x^{k} } \right] = 0, \end{gather}$$

(2.41c)$$\begin{gather}\zeta_0 (x) \frac{ \partial \hat{\vartheta} }{ \partial t } + \sum^{k=3}_{k=0} \left[ \eta_k (x) \frac{ \partial^{k} \hat{h} }{ \partial x^{k} } + \mu_k (x) \frac{ \partial^{k} \hat{q} }{ \partial x^{k} } + \nu_k (x) \frac{ \partial^{k} \hat{\vartheta} }{ \partial x^{k} } \right] = 0, \end{gather}$$

where $( \hat {q} , \hat {h} , \hat {\vartheta } )$ are infinitesimal disturbances to the steady-state solutions $( q_s, h_s, \vartheta _s )$ and $\varphi _k = \{ \alpha _k , \beta _k , \gamma _k , \xi _k , \zeta _k , \eta _k , \mu _k , \nu _k \}$ are linearised periodic coefficients – the latter can be found in the online supplementary material.

It is important to note that Squire's theorem (Squire's theorem states that two-dimensional instabilities are the least stable modes) is not always valid in the presence of thermo-capillarity (Kalliadasis et al. Reference Kalliadasis, Demekhin, Ruyer-Quil and Velarde2003a); furthermore, long waves are also not necessarily the most unstable in film flow down wavy substrate (Pollak & Aksel Reference Pollak and Aksel2013). This means three-dimensional disturbances of all possible wavenumbers should be considered. However, since the present formulation is strictly only applicable to small $Ma$, when Squire's theorem still applies, only two-dimensional instabilities are considered here. The disturbances are modelled by a sum of Floquet wave harmonics, which take the form

(2.42a)$$\begin{gather} \hat{q} = \textrm{e}^{- \textrm{i} \omega t} \,\textrm{e}^{2 {\rm \pi}\textrm{i} Q x} \sum^{m=\mathbb{F}}_{m={-}\mathbb{F}} \bar{q}_m \,\textrm{e}^{2 {\rm \pi}\textrm{i} m x}, \end{gather}$$

(2.42b)$$\begin{gather}\hat{h} = \textrm{e}^{- \textrm{i} \omega t} \,\textrm{e}^{2 {\rm \pi}\textrm{i} Q x} \sum^{m=\mathbb{F}}_{m={-}\mathbb{F}} \bar{h}_m \,\textrm{e}^{2 {\rm \pi}\textrm{i} m x}, \end{gather}$$

(2.42c)$$\begin{gather}\hat{\vartheta} = \textrm{e}^{- \textrm{i} \omega t} \,\textrm{e}^{2 {\rm \pi}\textrm{i} Q x} \sum^{m=\mathbb{F}}_{m={-}\mathbb{F}} \bar{\vartheta}_m \,\textrm{e}^{2 {\rm \pi}\textrm{i} m x}, \end{gather}$$

where $\omega$ is the angular frequency, $Q \in [ 0 , 1]$ is the Floquet parameter (i.e. wavenumber), $\mathbb {F} \in \mathbb {Z}$ is the number of Floquet harmonics and $( \bar {q}_m , \bar {h}_m , \bar {\vartheta }_m )$ are the disturbance amplitudes. The periodic coefficients $\varphi _k (x)$ depend upon the steady-state solutions which must be found numerically prior to a stability analysis. The exception is the steady-state flow rate $( q_s )$ which according to (2.26) will be constant. A suitable approximation for $q_s$ is found by assuming the liquid film is thick enough for the flow rate to remain close to the flat-film solution, i.e. the Nusselt solution (Scholle, Wierschem & Aksel Reference Scholle, Wierschem and Aksel2004); putting $h = 1$ in (2.35a) yields $q_s = 2 / 3$.

Substituting (2.42) into (2.41) and applying a Fourier transform

(2.43)\begin{equation} \hat{\boldsymbol{f}} (n) = \int^{1}_{0} \boldsymbol{f} (x) \,\textrm{e}^{{-}2 {\rm \pi}\textrm{i} n x} \,\textrm{d}\kern0.7pt x \quad \text{for} \ n ={-} \mathbb{F} ,\ldots, \mathbb{F} , \end{equation}

leads to the generalised eigenvalue problem $( \boldsymbol {A} - c \boldsymbol {B}) \boldsymbol {\hat {x}} = 0$ involving the phase velocity $c = \omega / ( 2 {\rm \pi}Q )$ and the eigenvector $\boldsymbol {\hat {x}} = ( \bar {q}_n , \bar {h}_n , \bar {\vartheta }_n )^{\textrm {T}}$.

Focusing on the temporal stability, $Q \in \mathbb {R}$ with $c \in \mathbb {C}$; the problem then reads

(2.44a)\begin{equation} \sum^{m=F}_{m={-}F} \delta_{m, n} \left[ - c_{} \bar{h}_m + \left( 1 + \frac{m}{Q} \right) \bar{q}_m \right] = 0, \end{equation}

(2.44b)\begin{align} &\sum^{m=F}_{m={-}F} \left[ - c ( 2 {\rm \pi}\textrm{i} Q ) \sum^{k=2}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} \hat{\alpha}_{k, n - m} \bar{q}_m + \sum^{k=2}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} \hat{\xi}_{k, n - m} \bar{\vartheta}_m \right. \nonumber\\ &\quad \left.+ \sum^{k=4}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} [ \hat{\beta}_{k, n - m} \bar{h}_m + \hat{\gamma}_{k, n - m} \bar{q}_m ] \right] = 0, \end{align}

(2.44c)\begin{align} &\sum^{m=F}_{m={-}F} \left[ - c ( 2 {\rm \pi}\textrm{i} Q ) \hat{\zeta}_{0, n - m} \bar{\vartheta}_m + \sum^{k=3}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} \hat{\eta}_{k, n - m} \bar{h}\right. \nonumber\\ &\quad \left.+ \sum^{k=1}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} \hat{\mu}_{k, n - m} \bar{q}_m + \sum^{k=2}_{k=0} (2 {\rm \pi}\textrm{i})^{k} (Q + m)^{k} \hat{\nu}_{k, n - m} \bar{\vartheta}_m \right] = 0, \end{align}

for $n = - F ,\ldots, F$; where $\delta _{m, n}$ is the Dirac delta function, $\hat {\varphi }_{k, n - m} = \int ^{1}_{0} \varphi _k (x) \,\textrm {e}^{- 2 {\rm \pi}\textrm {i} ( n - m) x} {\textrm {d} x}$ are the Fourier expansion coefficients and the entries of the matrices $\boldsymbol {A}_{m, n}$ and $\boldsymbol { B}_{m, n}$ are constructed using the periodic coefficients $\varphi _k (x)$. In the work reported here, the eigenvalues $c$ were found using Matlab's built-in subroutine eig and the stability determined by the eigenvalue with the largest positive imaginary part, i.e. the largest growth rate; if there are no eigenvalues with a positive imaginary part then the film is considered stable. Neutral stability is defined as when instabilities neither grow nor decay in an exponential fashion (the centre manifold of the system); this ensues when the eigenvalue of the most unstable mode is wholly real.

When performing the analysis, it is sufficient to consider only half the interval of the Floquet parameter, $Q \in [ 0 , \frac {1}{2} ]$. This is because of the symmetry, $c_n ( - Q ) = c^*_n ( Q )$, and periodicity, $c_n ( Q + 1 ) = c_n ( Q )$, of the eigenvalues; hence, $c_n ( \frac {1}{2} + Q ) = c^*_n ( \frac {1}{2} - Q )$.