Robust low-dimensional modelling of falling liquid films subject to variable wall heating

Accurate low-dimensional models for the dynamics of falling liquid films subject to localized or time-varying heating are essential for applications that involve patterning or control. However, existing modelling methodologies either fail to respect fundamental thermodynamic properties or else do not accurately capture the effects of advection and diffusion on the temperature profile. We argue that the best-performing long-wave models are those that give the surface temperature implicitly as the solution of an evolution equation in which the wall temperature alone (and none of its derivatives) appears as a source term. We show that, for both flat and non-uniform films, such a model can be rationally derived by expanding the temperature field about its free-surface values. We test this model in linear and nonlinear regimes, and show that its predictions are in remarkable quantitative agreement with full Navier–Stokes calculations regarding the surface temperature, the internal temperature field and the surface displacement that would result from temperature-induced Marangoni stresses.

FIGURE 1. (Colour online) Concept for thermocapillary control system. Near-real-time observations of the surface deformation are used to activate localized wall heating strips, in such a way that the system is driven to stability or instability.
system. However, practical implementation of feedback via blowing and suction through the supporting wall presents significant challenges in terms of delivering well-controlled, real-time inputs, and would likely require invasive modifications of any industrial process. A more appealing method of applying feedback is presented by the possibility of applying time-dependent and spatially varying heating of the solid wall, again in response to observations, using a system such as illustrated in figure 1. In principle, implementation could require only a small modification of the heating elements currently used in experiments, replacing a uniformly heated wall with strips that are controlled electronically and switched frequently between different states. The applied heating would affect the hydrodynamic flow through the temperature dependence of physical properties such as surface tension.
Owing to its appearance in the mass conservation equation, blowing and suction of fluid has a straightforward effect on film dynamics (Thompson, Tseluiko & Papageorgiou 2016b), leading to the success of a simple control scheme where fluid is injected (resp. withdrawn) in regions where the film is thin (resp. thick). Interaction with the background flow and film non-uniformities lead to relatively small modifications of this basic picture, and so the success of control is largely model-invariant. However, heating applied via localized elements embedded in the wall has a much more subtle effect on film dynamics: due to advection and diffusion, we expect that temporally transient heating applied via a localized heating strip will influence the fluid temperature mainly downstream, over a wider region than the original strip, and persisting over some time. We can quantify the importance of advection relative to diffusion via the Péclet number Pe = h s U s /κ, where h s is a typical film thickness, U s is a typical free-surface velocity and κ is the thermal diffusivity. The ratio between Pe and the Reynolds number Re = h s U s /ν, where ν is the kinematic viscosity, is the Prandtl number, Pr = ν/κ = Pe/Re, which is a material property of the fluid. For water, Pr ≈ 7, meaning that advection of heat is significant whenever inertia is, so non-negligible for films with inertia-driven instabilities.
As described below, there has been considerable attention paid to modelling and simulating wave dynamics in the case of uniform wall heating. However, there has been relatively little focus on the specific challenges that heterogeneous heating entails. Our aim in this study is precisely to systematically develop robust long-wave models to quantitatively predict the spatial and temporal distribution of temperature 846 A. B. Thompson, S. N. Gomes, F. Denner, M. C. Dallaston and S. Kalliadasis which arises in response to particular specified wall inputs. We will focus on an idealized system in which we are able to specify the temperature at the wall, Θ(x, t), exactly. For most fluids, viscosity is strongly dependent on temperature, but, like in nearly all other studies reviewed below, we will follow the standard assumption of temperature-independent viscosity. We will therefore assume that the fluid temperature only affects the hydrodynamic flow through the variation of surface tension with temperature, so that the only coupling of heat upon the flow is via thermocapillarity and Marangoni stresses. As a result, we seek to model the surface temperature S(x, t) and its dependence on wall temperature. We are thus able to direct our attention to the effects of advection and diffusion in shifting, delaying and smoothing the inputs that we apply.
For isothermal falling liquid films subject to inertia, there have been a number of different long-wave modelling approaches, including those by Benney (1966), Shkadov (1969) and Ruyer-Quil & Manneville (2000). The Benney model involves a single evolution equation for the film height h(x, t). This model is easy to use and accurately captures the transition to instability, but displays unphysical blow-up behaviour (Pumir, Manneville & Pomeau 1983). Shkadov (1969) developed an integral-boundary-layer model, involving coupled evolution equations for the film height and the mean velocity. While blow-up is avoided, it does not predict well the instability threshold compared to the full system except for large inclination angles (Benjamin 1957;Yih 1963). The weighted-residual methodology developed by Ruyer-Quil & Manneville (2000) combines the best features of both models: two or more coupled evolution equations are used (describing film height and downstream flux, among others), the critical Reynolds number is accurately captured, and blow-up does not arise. Furthermore, when compared to experiments and numerical solutions of the Navier-Stokes equations, the reduced weighted-residual model with second-order viscous but not inertial terms accurately captures the maximum and minimum height of solitary waves even up to Re as large as 115 in our scaling (Denner et al. 2018).
Previous research on non-isothermal films has focused on two fundamental problems, exploring firstly the effect of a uniformly heated bounding wall on the linear stability and nonlinear wave development of the film, and secondly the effect of imposed steady, spatially localized heating on the surface shape and stability. Marangoni stresses only arise if the surface temperature is non-uniform. For a uniformly heated wall, such surface temperature variation requires non-uniform film thickness, as is relevant to travelling waves and other time-dependent behaviour. The leading-order long-wave model to relate the surface temperature S(x, t) to the wall temperature Θ(x, t) is S(x, t) = Θ 0 /(1 + Bi h(x, t)), where the Biot number Bi 0 is used to parametrize heat loss at the free surface. The implied dependence of S on h via (1.1) can lead to a new mode of instability. If BiΘ 0 > 0, the surface is hottest when nearest to the heated wall and cooler when further away. The temperature affects the film evolution via Marangoni stresses, with a contribution to the mass flux in the direction from hot to cold. In combination, provided that BiΘ 0 > 0, a heated wall is destabilizing, while a cooled wall is stabilizing.
The leading-order thermal model (1.1) completely omits all advection and axial diffusion effects. In fact, the axial diffusion of heat is significantly enhanced by advection: this is the Taylor-Aris dispersion effect, whereby advection and cross-stream diffusion combine to yield a large effective streamwise diffusion (Taylor 1953;Aris 1956). Advection and streamwise diffusion enter long-wave models as first-and second-order corrections, respectively. Within the context of a uniformly heated wall and a non-uniform film, relevant models include those by Ruyer-Quil et al. (2005) in a gradient expansion, while Trevelyan & Kalliadasis (2004), Trevelyan et al. (2007) and Thiele, Goyeau & Velarde (2009) used a weighted-residual methodology. In their book, Kalliadasis et al. (2012) derive various second-order models for heated falling films, including Marangoni and inertial corrections to the parabolic velocity field, and note that there is a whole family of long-wave models that are equivalent at up to second order, with no obvious criteria to distinguish amongst them.
Long-wave models for heat transfer do not necessarily predict temperature distributions that are in accordance with the laws of thermodynamics. In particular, it is imperative that the fluid temperature should remain bounded between the extremal values of the wall and air temperature, and that there should be no local extrema of temperature inside the fluid layer. However, all but the simplest long-wave models predict negative (relative to both ambient and wall) temperatures in at least some parameter regimes. Attempting to remedy this issue, Chhay et al. (2017) recently developed a new model for heat transfer from a uniformly heated, vertical wall, for which viscosity and surface temperature are assumed to vary with temperature. Their model is based on Saint-Venant-type averaging, and comprises a first-order hyperbolic equation for the surface temperature, coupled to evolution equations for film height and the depth-averaged velocity. This new model delays the onset of unphysical negative temperatures at small Péclet numbers; the authors speculate that it would be impossible to develop an accurate long-wave model for non-uniform films that never predicts a negative temperature. Although our analysis here does not focus on uniformly heated walls, thermodynamic consistency is one of our main criteria for a physically reasonable long-wave model.
For films subject to non-uniform wall heating, the canonical problem concerns the effect of a localized heat source, i.e. a single hot wire or heating strip, embedded in the spanwise direction in an otherwise insulating or fluid-cooled rigid wall. A pronounced capillary ridge forms above the heated region, which becomes unstable to spanwise perturbations when the Marangoni number exceeds a critical value; the instability manifests as steady rivulets, equally spaced in the spanwise direction. Kabov, Marchuk & Chupin (1996) used detailed measurements of surface temperature and height profiles, combined with solution of the Nusselt problem and a heat transfer problem within the substrate, to infer the imposed heat flux at the free surface. Rivulets are associated with sharp gradients in surface temperature, of the order of 10 K mm −1 . The combination of large Prandtl number and moderate Reynolds number suggests that advection should be significant, supported by the slow downstream decay of the observed surface temperature. Surface velocity measurements performed using tracer particles show that the ridge instability is associated with the development of a stagnation point on the surface, corresponding to a backflow region within the fluid (Kabov et al. 2001(Kabov et al. , 2002. Marchuk & Kabov (1998) simulated the two-dimensional (2-D) problem for ridge formation, using a lubrication equation at zero Reynolds number for the hydrodynamic problem subject to a known surface temperature, iteratively coupled to a 2-D numerical solution of the heat equation using the flow field from the lubrication equation in order to determine the surface temperature. In this approach, inertial effects are neglected but advection is included, a combination suitable only for fluids with relatively large Prandtl number. The boundary conditions for the heat problem are specified flux within the heater, and either zero flux or specified temperature outside; the air-fluid surface is assumed to be perfectly insulating (Bi = 0). Comparison with the surface temperature and displacement measurements by Kabov et al. (1996) show satisfactory agreement. 848 A. B. Thompson, S. N. Gomes, F. Denner, M. C. Dallaston and S. Kalliadasis Frank (2003) performed three-dimensional (3-D) Navier-Stokes simulations using the wall temperature profile inferred by Kabov et al. (1996). The calculations, using calibrated temperature-dependent profiles for both fluid viscosity and surface tension, were in good agreement for the critical heat flux and spanwise wavelength for the onset of rivulets. Frank (2003) also found that spanwise instability arises even if Bi = 0, indicating that a different mechanism is at work than in the Marangoni-induced instability in a uniformly heated system. Frank & Kabov (2006) conducted experiments and simulations for films across a wider range of parameters. They showed numerically that the decay rate of surface temperature downstream is sensitive to the Biot number, with a comparison to experimental measurements of this decay rate suggesting that Bi 0.03. Comparisons between 3-D numerics and experiments for the critical heat flux and corresponding spanwise wavelength at the onset of rivulets are in good agreement.
There have been several attempts to use a long-wave framework to replicate the experimental results for the development of a capillary ridge, its spanwise stability, and the selection of nonlinear patterns after the onset of instability. In order to adequately describe this problem, a long-wave model is required to relate the surface temperature to the wall temperature as mediated by a flow that is in principle neither uniform nor steady. Kalliadasis, Kiyashko & Demekhin (2003) examined the response of a falling liquid film to a single wide heating strip. They derived both the hydrodynamic and energy equations using an integral-boundary-layer approximation, eventually resulting in an evolution equation for the surface temperature involving up to first derivatives with respect to space and time of the wall temperature and surface temperature. This model is accurate to first order, capturing advection and not requiring a flat film, but does not include any axial diffusion terms. The model indeed predicts that the film surface exhibits a bump in regions of positive gradient of the wall temperature, and that this bump becomes unstable to transverse perturbations when the Marangoni number is sufficiently large. However, no quantitative comparison with experiments is attempted.
At around the same time, Skotheim, Thiele & Scheid (2003) used long-wave modelling and experiments to investigate the capillary-ridge problem. They completely bypassed the thermal problem by specifying the surface temperature a priori, while a very low Reynolds number means that inertial effects are also neglected. The ridge shape and wavenumber at onset of instability are compared to experiments, with the Marangoni number being treated as a fitting parameter in this comparison, chosen to ensure the same bump height in experiments and theory. The results are in reasonable agreement for the ridge shape and spanwise rivulet wavenumber, but the critical Marangoni number for the onset of rivulets is around five times too high. It is noteworthy that in their study the stability analysis is performed under the assumption of a periodic array of heaters; results for varying periodicity are extrapolated to infer the effect of a single heater. Scheid et al. (2002) explored the effect of imposed steady sinusoidal heating on travelling waves and showed that, for sufficiently large Marangoni number, the travelling waves become 'frozen', corresponding to a steady state induced by the patterning. They worked at moderate Reynolds number, using a Benney equation for the hydrodynamic problem and, like Skotheim et al. (2003), only the leading-order model for the thermal problem (this combination is suitable only for fluids with small Prandtl number); they also briefly discuss the effect of imposing a heat flux at the wall or a mixed heat boundary condition at the wall. To mimic the expected effects of advection and diffusion, a wall temperature profile is chosen with a steep ramp up to a uniform temperature followed by a slow downstream decay. Like Skotheim et al. (2003), the Marangoni number is chosen to ensure the same bump thickness at the onset of instability. The results are in broad agreement with the experimental results from Scheid et al. (2000), with discrepancies attributed to the unknown Biot number and the lack of temperature-dependent viscosity in the model. D'Alessio et al. (2010) considered uniform heating on a wavy wall. The focus of their study is on non-uniform films with non-trivial bottom topography, but the heat equation used is more sophisticated than in many studies. The starting point is an assumed temperature profile, linear in the cross-film coordinate y. Such a profile cannot satisfy the heat flux condition on the free surface, but Biot-number effects are included via an integral condition, which also allows for the inclusion of an axial diffusion term. The resulting equation is a parabolic partial differential equation (PDE): an advection-diffusion equation for S(x, t) whose source term is Θ − (1 + Bi)S. We will show later that there is a unique model with this structure which is both accurate to second order in the long-wave limit when applied for flat films and also valid in the case of non-uniform heating. Comparing D'Alessio et al.'s (2010) equation in the case of flat topography and flat free surface with the equation we derive in § 4.4, we find that they have failed to capture the enhancement of diffusion by flow, and there are also some discrepancies in the coefficients.
Finally, Blyth & Bassom (2012) explored the use of non-uniform heating via the Marangoni effect to passively control the surface deformation resulting from a wavy wall. They worked in the limit of Re = Pe = 0, and solved the 2-D equations numerically without any long-wave approximations. For small topography amplitudes, they showed that sinusoidal heating can flatten the fluid surface, while large enough steady uniform heating induces the free surface to be in phase with the wall shape; ultimately, for a flat wall with a downward step, localized cooling could be used to eliminate a capillary ridge.
Our main focus in this paper is on formulating thermodynamically consistent long-wave models for heat transfer through a falling liquid film. The challenge is the accurate capture of advection and diffusion effects, which arise even in the case of a flat film subject to a parabolic velocity profile. We begin by stating the full 2-D governing equations for both flow and heat transfer in § 2. In § 3, we consider the special case of heat transfer through a flat film with a parabolic velocity profile and no Marangoni coupling. This problem can be approached by Fourier transforms but has no closed-form solution in terms of elementary functions at non-zero Pe. Instead, we analyse the equations for each Fourier mode and prove five thermodynamic properties (P1-P5) that are obeyed by the full 2-D system and will serve as qualitative tests for long-wave models.
In § 4, we derive four long-wave models, each accurate to second order in the long-wave limit for the case of a flat film subject to non-uniform heating. These are a Benney equation ( § 4.1), a weighted-residual equation ( § 4.2), a new model that we call 'projected' ( § 4.3) and a second-order parabolic model with similarities to an advection-diffusion equation ( § 4.4). Qualitative performance against the five properties is discussed, with the new projected model and the parabolic model the most promising. In § 5, we quantitatively compare the four long-wave models with results from the full 2-D system for sinusoidal and strip-like heating. We find that only the projected model gives good agreement across a wide range of conditions, and subsequently we focus on this model alone. This model is a second-order, linear, hyperbolic equation, giving S(x, t) implicitly as a function of Θ. It changes character from subcritical to supercritical at a finite value of Pe. In the remainder  of § 5, we explore the implications of this change in criticality for S(x, t) and for the corresponding internal temperature field. We also investigate the response to time-dependent heating. In § 6, we consider the extension of the new projected model to the case of a non-uniform film. Steady states for the film height and surface temperature at both small and large deformation are presented in § 7, with a comparison to computational fluid dynamics solutions of the full, nonlinear 2-D Navier-Stokes system. We draw our conclusions in § 8.

Dimensional equations
We consider a 2-D flow, governed by the incompressible Navier-Stokes equations and with all bulk properties of the fluid (dynamic viscosity µ, density ρ, thermal conductivity k and specific heat capacity C p ) independent of temperature. We use axes aligned with the wall, as illustrated in figure 2; the coordinate x and velocity component u are in the downslope direction, while the coordinate y and velocity component v are perpendicular to it. The free surface is defined by y = h(x, t), and the wall is inclined at an angle β to the horizontal. We will suppose that the air temperature remains constant at a temperature θ A , set θ A + θ (x, y, t) to be the temperature within the film, and θ A + Θ(x, t) to be the wall temperature.
Heat has no effect on the bulk properties, so within the fluid layer we have the usual momentum equations, Robust low-dimensional modelling of heated falling films 851 where p is the fluid pressure and g the acceleration due to gravity, along with conservation of mass in the form Within the fluid layer, the temperature θ(x, y, t) evolves according to the heat equation On the rigid wall at y = 0, the fluid velocity vanishes, and the fluid temperature θ is prescribed: In the bulk equations, the velocity field can influence the temperature evolution, but there is no reverse coupling. The only feedback from the fluid temperature to the hydrodynamic fields is through the dynamic boundary condition, where σ is the stress tensor, n and t are unit vectors normal (directed into the air layer) and tangential to the free surface, respectively, and γ is the temperaturedependent coefficient of surface tension. For simplicity, we will model this temperature dependence by a linear relation, as appropriate for small variations in temperature, but a nonlinear function could be used without difficulty. Hence we write where the positive constants γ A andγ , respectively, are the coefficient of surface tension of the fluid at the temperature of the air and the strength of surface tension variation.
The hydrodynamic system also satisfies the kinematic boundary condition, and we will assume a simple linear relation for the heat flux into the air, where the quantity L has dimensions of length. The relationship (2.9) is commonly used in the heat transfer literature, and follows from the assumption that the air layer is always held at θ = 0, combined with Newton's law of cooling.
2.2. Non-dimensionalization We will non-dimensionalize the system based on a typical film height h s and the corresponding surface speed U s = ρgh 2 s sin β/(2µ). Hence we define Thompson, S. N. Gomes, F. Denner, M. C. Dallaston and S. Kalliadasis We will scale θ, Θ and S on a typical temperature scale Θ, which corresponds to a characteristic temperature variation between the heated wall and the ambient air. The bulk equations become (dropping hats), each holding in the region 0 < y < h(x, t): and u x + v y = 0. (2.14) The dimensionless parameters in the bulk equations are the inclination angle β, and the Reynolds, Péclet and Prandtl numbers, defined here as The uniform state of an isothermal falling liquid film becomes unstable when Re > 1.25 cot β, and the initial instability is a long-wave variety. The Péclet number is related to the Reynolds number by the Prandtl number, which is a material constant, and is the ratio of diffusivities of momentum ν and heat κ = k/(ρC p ) in the fluid. For water, Pr ≈ 7. This means that, if the Reynolds number is large enough to induce instabilities, then advection cannot be neglected when considering heat transfer through the film, and so for consistency we should consider regimes with moderate Péclet numbers. In component form, the dynamic boundary condition becomes .16) and where Ca and Ma are the capillary and Marangoni numbers, respectively, The kinematic condition is unchanged in these new variables, and the heat loss condition becomes n · ∇θ = −Bi θ, on y = h(x, t), where Bi = h s /L 0 is the Biot number. The Biot number condition is a modelling convenience, bypassing the conjugate heat transfer problem with the air and, given the low thermal conductivity of air, is often taken to be small.

Fundamental properties of heat transfer in a flowing, uniform film
We begin our analysis by considering heat transfer in the full 2-D heat equation for the case where the film is flat and the flow is uniaxial: We note that, if Ma = 0, this flow field exactly satisfies the governing equations (2.11), (2.12) and (2.14) and the boundary conditions (2.17) and (2.16), regardless of θ . As the heat equation is linear in temperature, and the base flow is translationally invariant with respect to x and t, we can use a Fourier decomposition to analyse the temperature transmission, with the x and t dependence of all temperature fields taking the form exp(ikx + λt): with boundary conditions The general solution to (3.3) involves parabolic cylinder functions as a result of the quadratic terms in u(y). However, this solution is not particularly illuminating, and of course the boundary-value problem for T(y) can be solved numerically. For our purposes, the qualitative behaviour of the solutions is more useful. The heat transfer behaviour must obey three intuitive properties for individual Fourier modes: P1. For steady, non-uniform heating, the amplitude of heat variation at the film surface is not larger than the amplitude of heat variation at the wall.
Proof. Let T(y) = r(y) exp(iφ(y)) where both r 0 and φ are real. As the system is steady, we can decompose the conditions above into Hence r(1) r(0)/(1 + Bi), and as Bi 0, we have r(1) r(0) which is exactly the condition required.
P2. At high wavenumber, the amplitude of heat variation at the surface tends to zero.
P3. If the wall temperature is held fixed at Θ = 0, the fluid temperature will always evolve towards θ = 0, regardless of the initial temperature distribution.
Proof. We write the eigenvalue λ = λ R + iλ I . Using the polar decomposition of T(y) as above, the system becomes (3.9) Then using the integral (3.6), As r(y) 0 (and is not zero everywhere), we conclude that Peλ R + k 2 0, and so λ R 0 for Pe > 0.
In terms of physical interpretations, P1 follows from the thermodynamic principle that heat must flow from hot to cold regions, so the heat flux should not lead to unforced hot spots. P2 relates to the dominance of diffusion at high k; in the absence of flow, axial diffusion means that r(1)/r(0) decays exponentially with k, and we would expect this decay rate to be enhanced by the stretching effect of the fluid flow, which leads to Taylor-Aris dispersion. P3 means that the constant temperature distribution (equal to the wall and air temperature) is stable.
The three properties above are clearly not exhaustive, but they are some of the most important. These three properties are derived via Fourier transforms under the assumption that the film surface is flat and the velocity field is purely Nusselt flow. Hence they only apply to uniform films, and require that Ma = 0.
We can also appeal to two properties that should hold for steady films, which could be uniform or non-uniform. For the sake of argument, we envisage periodic solutions, so that the maximum and minimum of the wall and surface temperatures are attained, and we do not have to consider behaviour as x → ±∞. These properties are as follows: P4. In a steady state, there are no local maxima or minima of the temperature inside the fluid layer.
Proof. This property is a special case of Hopf's maximum principle for linear elliptic PDEs (see e.g. Evans 2010). In steady state, the temperature equation is (3.12) At a local extremum of the temperature, ∇θ = 0, so the left-hand side of (3.12) is zero, and hence ∇ 2 θ = 0 at that point, so a simple maximum or minimum is not permitted. At the next order, taking the Laplacian of (3.12) yields If ∇∇ 2 θ = 0, then the stationary point is a saddle point, regardless of the value of ∇ 4 θ. Alternatively, if ∇∇ 2 θ = 0 at the point in question, then ∇ 4 θ = 0 also. The process of taking the Laplacian can be repeated as needed, but the left-hand side will always contain odd derivatives of lower order than the right-hand side. This implies that the first non-zero derivative cannot be of even order, and hence that the temperature field cannot have an internal local maximum or minimum.
P5. The surface temperature is bounded by the maximum and minimum of the wall temperature and the air temperature.
Proof. We start with the case Bi > 0. The boundary condition at the surface (2.19) gives that n · ∇θ = −Bi S, where n is an outward-pointing unit normal. Consider the point on the surface with maximum value S max . If S max > 0, then the boundary condition (2.19) implies that the outward normal derivative of θ is negative. Hence there must be a point x * , just below the surface, such that θ (x * ) > S max . Since θ at this interior point must be bounded by the maximum value on the boundary of the fluid region, and is greater than the maximum on the free surface, it must be bounded by the maximum value on the wall: thus (3.14) The corresponding argument for the minimum surface temperature will give S min min{0, Θ min }.
For the case Bi = 0 with a flat film, we can proceed by mirroring the fluid region in the free surface, so that points on the surface become interior points of an extended domain. For a flat surface of height H, we can extend u, v and θ into the region H y 2H according to Given that v(H, y) = 0 and θ y (H, y) = 0 from the boundary conditions, we find that u = (u, v) is continuous and θ is smooth over the extended domain. In this domain, θ satisfies the same PDE as above, with boundary conditions θ ( Thus all points of the original free surface y = h are interior points of the new problem, and must therefore be bounded between Θ min and Θ max .
We seek long-wave models accurate to second order in the long-wave limit, and consistent with as many of the properties listed above as possible. In practice, we find that there is a unique model (in the sense discussed in § 4.4) that always obeys P1-P3, but even for this model it is possible to construct scenarios where P4 and P5 are not satisfied.

Long-wave models for heterogeneous heat transfer through a uniform Nusselt flow
In this section, we derive a number of long-wave models, each a 'local' PDE for the surface temperature, requiring only one spatial dimension and equivalent at up to second order in the long-wave regime. The models have qualitative differences in their structure (explicit, hyperbolic, parabolic), which means that they perform very differently when compared to the qualitative properties (P1-P5) listed above.
We will require all models to be consistent to second order in the long-wave limit with the full 2-D system for heat transfer in a flat film. In this long-wave derivation, we set X = δx and T = δt, with δ 1 and X, T = O(1). The heat equation is linear, so we do not need to set a scale for θ(X, y, T). In these new long-wave variables, we obtain and θ y + Bi θ = 0, on y = 1. (4.3) Our aim is to retrieve S(X, T) = θ(X, y = 1, T).

4.1.
Regular expansion in long-wave limit: a Benney-type model The most obvious long-wave model is to obtain a small-wavenumber expansion about the known linear temperature profile for zero wavenumber (k = 0 or δ = 0). The resulting expansion is regular, and is derived in much the same way as the Benney model for thin-film hydrodynamics. We pose an expansion of θ (X, y, T) for small δ: (4.4) The PDE (4.1) yields the sequence of problems: The boundary conditions come from (4.2) and (4.3) with the wall temperature Θ(X, T) treated as a known O(1) quantity, so that with θ ny = −Bi θ n at y = 0 for n = 0, 1, 2, . . . . We find that θ 0 , θ 1 and θ 2 are respectively first-, fifth-and ninth-order polynomials in y. With θ 0 , θ 1 and θ 2 known, the surface temperature is simply S = θ(y = 1).
Retaining terms in the surface temperature up to and including O(δ 2 ), and then returning to the original variables, this model gives Thus we obtain S explicitly through knowledge of Θ and its derivatives with respect to space and time. We can write (4.7) schematically as where L B is a constant-coefficient, second-order linear operator. If Θ = 0, then S = 0, so property P3 is satisfied: in fact, the temperature distribution tends to zero instantaneously, which is not physically realistic. However, a more important problem arises for steady transmission: (4.7) is quantitatively accurate for small k, but diverges for large k, where the amplitude grows as k 2 rather than decaying exponentially. If we were to calculate further corrections in the long-wave limit to obtain an nth-order-accurate model, the amplitude would grow as k n , and hence further long-wave corrections would lead to increasingly rapid divergence when evaluated in the short-wave limit. This means that the model fails properties P1 and P2. This explicit model is therefore ill-posed at short wavelengths, and of little use for finding S(x, t) in practical cases. However, the explicit expansion (4.7) is useful for testing quantitative agreement in the strictly long-wave limit, and might also have some relevance for inverse problems, where we wish to determine Θ(x, t) to deliver a particular surface temperature.

4.2.
Weighted-residual approach for temperature equation Inspired by the success of the weighted-residual approach developed by Ruyer-Quil & Manneville (2000) in modelling thin-film hydrodynamic flow at moderate Re by coupled evolution equations for h(x, t) and q(x, t), a number of authors have similarly developed weighted-residual models for heat transfer through a flowing film. Previous papers usually present the models correct only to first order in the long-wave parameter δ (but note that, in their appendix, Kalliadasis et al. (2012) list a full second-order model, including Marangoni and inertial corrections to the velocity field). Furthermore, nearly all studies assume that heating is uniform or steady.
In a typical derivation of a weighted-residual temperature equation, the temperature is written in terms of basis functions, e.g. . (4.9) Note that, in contrast to the hydrodynamic case, the individual functions ψ n do not satisfy the heat loss boundary condition at y = 1. The unknown coefficients a n are determined by requiring that the temperature must satisfy a combination of strong and weak conditions. Strong conditions always include that the temperature or flux matches its prescribed value at the wall. Typical weak conditions are derived by 858 A. B. Thompson, S. N. Gomes, F. Denner, M. C. Dallaston and S. Kalliadasis multiplying the heat equation by a suitable weighting function and integrating across the depth of the film. As noted by Trevelyan et al. (2007) and Kalliadasis et al. (2012), even for a given order of long-wave accuracy, different choices of residual weights will lead to different model equations with different coefficients; the models will only necessarily be equivalent in the long-wave limit. In this section, we follow earlier studies by choosing weight functions for the nth-order weak conditions as y n−1 , but we apply only the first two weak conditions. We supplement the boundary conditions on θ at y = 0 and y = 1, with two strong conditions: we require the residual R(X, y, T), defined as to be equal to zero at both the wall y = 0 and the fluid surface y = 1. With the two weak conditions, the bulk residual yields four equations in total: We expand the temperature field as a polynomial in y, with coefficients that depend on X and T: θ (X, y, T) = Θ(X, T) + a 1 (X, T)y + a 2 (X, T)y 2 + a 3 (X, T)y 3 + a 4 (X, T)y 4 + a 5 (X, T)y 5 .
(4.13) This expansion automatically satisfies θ = Θ at y = 0; we wish to obtain an equation for S(X, T) = θ(X, 1, T). Five unknown coefficients are needed to obtain a model accurate to second order in the long-wave limit. We start by applying (4.10), which allows us to write a 1 and a 2 in terms of S, a 3 , a 4 and a 5 ; likewise (4.12a) does not require any derivatives of the unknown coefficients, allowing the algebraic elimination of a 3 . In contrast, the final three equations, (4.12b-d), each involve derivatives of S(x, t), a 4 (x, t) and a 5 (x, t). These three conditions yield coupled linear evolution equations for S, a 4 and a 5 , each first-order in time. However, as the linear operators commute, we can take a suitable linear combination to obtain an equation linking S and Θ alone. This equation features up to third-order derivatives in time, and sixthorder in space, of both S and Θ.
Once we have obtained a single evolution equation for S(x, t), we discard all terms in this third-order evolution equation that are smaller than O(δ 2 ) (other terms at O(δ 3 ) having already been neglected in our analysis), yielding the following governing equation: ( (4.14) The weighted-residual equation (4.14) involves second derivatives with respect to time and space of both S and Θ; its compact form is where both L W1 and L W2 are second-order linear operators. For steady heating, the solution to (4.14) is in agreement with (4.7) for small δ at up to second order. However, solutions to (4.14) do not always obey property P1 (surface amplitude always smaller than wall amplitude) and may also violate P2 (surface amplitude vanishes at short wavelength). For example, when Pe = 6 and Bi = 0, the coefficient of S xx in (4.14) is close to zero and S ∼ −1.63Θ as k → ∞, thus violating properties P1 and P2. For general parameter values, the steady heat transmission amplitude tends to a non-zero constant for large k (governed by the coefficients of the Θ xx and S xx terms), again contradicting property P2. However, all temporal eigenvalues λ have negative real part, and so the state S = 0 is stable if Θ = 0; hence P3 is satisfied.

Projection approach
The Benney model yields S explicitly as a function of Θ and its derivatives. In contrast, the weighted-residual model gives S implicitly by solving an evolution equation, with a source term that involves Θ and its derivatives. Both models are equivalent in the long-wave limit, but the weighted-residual model is better behaved in the short-wave limit in that the surface temperature is bounded. Here we go one step further, seeking an evolution equation for S in which only Θ, and not its derivatives, is the source term, so that the schematic form should be where L P is a second-order linear operator. Such a model can be derived as the 'inverse' counterpart to the Benney equation. We pose an expansion for small δ: θ (X, y, T) = θ 0 (X, y, T) + δθ 1 (X, y, T) + δ 2 θ 2 (X, y, T) + · · · . (4.17) However, this time we treat S(x, t) as a known, O(1) quantity, determined by the longwave PDE (4.1) and the boundary conditions θ(X, 1, T) = S(X, T), θ y (X, 1, T) = −Bi S(X, T). (4.18a,b) We expand the PDE (4.1) and the two boundary conditions (4.18) for small δ to obtain the sequence of PDEs (4.5) with boundary conditions θ 0 = S, θ 0y = −Bi S at y = 1, and θ 1,2 = 0, θ 1,2y = 0 at y = 1. Given that u(y) = y(2 − y), the solutions for θ 0 , θ 1 and θ 2 are polynomials in y, with θ 0 involving S(X, T), θ 1 featuring S X and S T , and θ 2 involving first and second derivatives of S. With θ (x, y, t) known by truncating the small-δ asymptotic sequence after O(δ) 2 , we finally impose the boundary condition at the wall: θ = Θ(x, t) at y = 0 (i.e. projecting the expansion onto y = 0). This statement requires (1 + Bi)S + Pe 1 2 + For steady heating, it is possible to show that the model obeys both property P1 (surface amplitude always smaller than wall amplitude) and property P2 (surface amplitude vanishes at short wavelength). Furthermore, the amplitude of surface variation decays monotonically as the wavenumber increases. The temporal eigenvalues always have negative real part, and so the state S = 0 is stable, thus satisfying property P3 (if the wall is held at fixed temperature, the fluid temperature will always evolve towards zero). We will discuss properties P4 and P5 later (see § 5) as these do not relate to individual Fourier modes.

Other second-order models with special properties
We have highlighted above three long-wave models that each agree to second order in the long-wave limit. The Benney model and our new projected model are special in that they give S(x, t) and Θ(x, t) explicitly, respectively. The weighted-residual model is between these two extremes, and is one of a family of linear second-order models, of which the general form is aS + bS x + cS xx + dS t + eS xt + f S tt = Θ + BΘ x + CΘ xx + DΘ t + EΘ xt + FΘ tt . (4.20) This equation has 11 unknown coefficients, which should be constant in space and time, but could depend on Pe and Bi.
We can analyse the solution behaviour for slowly varying heating by considering solutions of the form with δ 1. The solution of (4.20) forŜ iŝ Expanding (4.22) for small δ, we obtain For (4.20) to be second-order-accurate for long-wave heating, we require the expansion (4.23) to agree with that from the Benney equation for terms up to δ 2 , for arbitrary K and Λ. Agreement only occurs if the coefficients match exactly, yielding six algebraic equations involving the 11 coefficients a, b, . . . , F, but up to five parameters that can be freely chosen. The six conditions for second-order consistency cannot be satisfied by models with only first-order derivatives; at least one second-order derivative is required. A particularly appealing model structure would be a diffusion-like equation, involving an S xx term but not S xt or S tt , and no second derivatives of Θ. Such a model can be derived in two ways: either require that D = E = F = e = f = 0 above, and choose the remaining coefficients to be second-order-accurate in the long-wave limit (this yields a unique solution for the coefficients); or differentiate the first-order counterpart of (4.19) with respect to x and t to eliminate S xt and S tt in the projected model. Either derivation yields the same result. For the particular case Bi = 0, the resulting equation is This is a second-order, parabolic evolution equation for S, asymptotically correct to second order, for both steady and time-dependent heating. When extended to non-zero Bi, the coefficients of (4.25) become rational functions of Pe and Bi (note that for no choice of Bi and Pe are any of the coefficients equal to zero). The parabolic model (4.25) satisfies properties P1-P3 for arbitrary Bi 0 and Pe 0. For Θ = e ikx , S/Θ decays as O(k −1 ) for long-wave heating, which is better than the weighted-residual model but worse than the new projected model. We will show later that this model does not obey property P5; we can construct examples with negative surface temperature, see figure 5.

Quantitative comparison of models
We have derived four long-wave models to describe the surface temperature of a film subject to controlled heating at the wall. The models agree to O(δ 2 ) when the solutions for S are expanded in the long-wave limit, but differ qualitatively in their predictions of behaviour in the short-wave limit. In this section, we consider quantitative comparisons of their predictions, for Fourier modes and also for a more physically realistic heating profile in the form of a smoothed heated strip. Again, we set Ma = 0, so that the hydrodynamic solution is simply Nusselt flow.

Fourier signature
The first test concerns the relationship between surface temperature and wall temperature for steady, sinusoidal heating. In this case, we set Θ = exp(ikx) and S =Ŝ exp(ikx). Figure 3 shows results forŜ as a function of k for the full 2-D system (3.3) and (3.4), the Benney model (4.7), the weighted-residual model (4.14), our new projected model (4.19) and the parabolic model (4.25). As expected, all results agree in the long-wave limit, but the Benney model rapidly diverges from the full for the long-wave models (4.7), (4.14), (4.19) and (4.25) and for the full 2-D system solved by finite-differencing discretization of (3.3) and (3.4). The Benney results are in poor agreement except at very small k. For Pe = 1, the other three long-wave models are all in good agreement with the full system, though we can see some deviation at large k for both the weighted-residual and parabolic models. At Pe = 10, these deviations have become more pronounced, but the new projected model is in remarkably good agreement with the full system across the full range of k.
2-D results as k increases. For Pe = 1, the projected, weighted-residual and parabolic models give similar results and are close to the full 2-D results, though at large k, the projected model is in best agreement. At Pe = 10, only the projected model is still close to the full 2-D results and is in fact in very good agreement over the whole range of k. We can also analyse the time decay of a fluid layer with a sinusoidal temperature distribution and an unheated wall, corresponding to the relaxation to a uniform state. Hence we set Θ = 0 and S =Ŝ exp(ikx + λt), and solve for λ to give non-trivial solutions. Results for each system are shown for the case that Bi = 0 and Pe = 1, 10 in figure 4. The Benney model is not shown, as decay is instantaneous. For each long-wave model shown in figure 4, all eigenvalues have negative real part and the eigenvalue corresponding to the slowest-decaying mode are in broad agreement (within 15 %) of the full system results at small k. As k increases, the temporal decay for each eigenmode becomes more rapid in all models. The full system has an infinite sequence of decaying eigenmodes for each k, while the parabolic model has one, and the weighted-residual and parabolic equations have two eigenvalues for every k. When comparing the decay rate of the slowest-decaying mode, the parabolic model is in the best quantitative agreement with the full system over a wide range of k. The long-wave predictions for the second eigenvalue are quantitatively very poor, but in both weighted-residual and projected models these modes are predicted to decay rapidly, and the eigenvalues for the first and second mode for a given model do not cross for any Pe, Bi and k. (Note that the Benney model has no eigenvalues, while the full system has infinitely many; the two shown here are the slowest-decaying.) As for figure 3, we take Bi = 0, with (a) Pe = 1 and (b) Pe = 10. At k = 0, the exact eigenvalues are transcendental (λPe = −(π/2 + nπ) 2 ), and hence cannot be in perfect agreement with any of our long-wave models. However, the long-wave eigenvalues with real part closest to zero are all within 15 % of the exact value. The long-wave predictions for the second eigenvalue are in very poor agreement, but in all cases this eigenmode decays rapidly and should not be dynamically significant.

5.2.
Response to steady strip-like localized heating Although there are significant differences between the models at large k, it is not obvious how important this will be in practice. We now consider the surface temperature that results from a localized, strip-like heat source, with width 2X and centred at x = 0, represented as in a periodic domain with period L = 2π/m. The sharpness of the strip is controlled by σ . We take σ = 10 and X = 2, so that the wall temperature is almost constant in a region of length 4, and positive but close to zero outside this region. The surface temperature as predicted by the full 2-D system and each long-wave model is shown in figure 5. We find that the Benney model gives nonsensical results for both Pe = 1 and Pe = 10, predicting significant regions with negative temperature, and very large deviations from the true results. At Pe = 1, the other three long-wave models are all in good agreement with the full system. However, for Pe = 10, the parabolic and weighted-residual predictions exhibit extra turning points. The parabolic model predicts a region with negative surface temperature near the upstream end of the heating strip, thus violating property P5. By contrast, the projected model is in good agreement with the full results at both Pe = 1 and Pe = 10. In the remainder of this study, we will focus our analysis on the new projected model (4.19), which has performed well at every test so far.
The coefficients a, b, d, e and f are always positive for Pe > 0 and Bi 0. In contrast, c changes sign at a critical value of Pe, which we will denote as Pe * . This critical value is a function of Bi (plotted in figure 6) and increases monotonically with Bi, from Pe * = 4.592 at Bi = 0 to Pe * = 6.563 as Bi → ∞. (Note that the weighted-residual equation (4.14) also has a critical Pe * , with 6.297 < Pe * (Bi) < 11.277. The parabolic equation (4.25) has no such sensitivity.) Viewed as a time-dependent PDE, (5.3) is hyperbolic if Pe > 0 and Bi 0. For Pe < Pe * , its two families of characteristics travel in opposite directions, while for Pe > Pe * , both characteristics travel downstream. Hence if a localized disturbance is introduced, at small Pe the influence will be felt both upstream and downstream, whereas at large Pe, there is no upstream effect.
The steady influence of a localized heating strip has a similar dependence on Pe. As a result, (5.3) is diffusion-dominated for Pe < Pe * in the sense that its two real characteristics travel in opposite directions, so that a localized heat source has both upstream and downstream influence on the surface temperature. For Pe > Pe * , the system is advection-dominated; both characteristics travel downstream, and a localized heat source has no upstream influence. At large Pe and large Bi, a third region arises (bordered by the condition b 2 − 4ac = 0), where both characteristics travel downstream, but the solutions for r in (5.5) are complex. It is in this region only that our model can predict negative surface temperatures arising from positive wall heating, potentially violating P5. temperature occurs at x = 0, in the middle of the strip. As Pe increases, the location of this maximum temperature moves downstream, the maximum value decreases, and the spatial decay is generally slower; in all cases S(x) 0 everywhere. Figure 7 also shows full 2-D results computed from (4.2) and (2.13) using Fourier transforms in x with a finite-difference discretization in y. The agreement between the full system and the new projected model is generally good. The most significant discrepancies occur when Pe is close to Pe * where the model becomes singular in that the coefficient of S xx vanishes, resulting in the surface temperature profile predicted by the long-wave model displaying an unphysical rapid change in gradient near the two ends of the smoothed heating strip. Variation of parameters allows us to construct solutions to the steady-state equation: (5.5a,b) and hence to demonstrate the positivity of the solution. The complementary function involves the two solutions of the form S = exp(rx), where r is defined in (5.5b). When Pe < Pe * , we have a, b > 0 but c < 0; the roots for r are both real but of opposite sign, so we can take r 2 < 0 < r 1 and the appropriate boundary conditions are that S → 0 as x → ±∞. Supposing that S = 0 except inside a strip −z < x < z, we can write the solution for general Θ(x) as x > z. (5.6) The temperature profile decays exponentially with respect to x outside the localized heating region. We know that c(r 2 − r 1 ) > 0, so if Θ(x) 0, we must have S(x) 0 everywhere. This proves property P5 for the case Pe < Pe * . If Pe > Pe * , both r 1 and r 2 have negative real part. With boundary condition S → 0 as x → −∞, the solution is If r 1 and r 2 are real, and Θ(x) 0 for all x, then the integral solution (5.7) means that S(x) 0 for all x. Hence for parameter values where r 1 and r 2 are real, regardless of Pe, we have shown that for arbitrary positive substrate heating, the surface temperature is also positive; thus addressing property P5. The only condition under which the projected model can predict negative surface temperature with non-negative wall temperature is if the roots for r are complex, as then the integrand in (5.7) is oscillatory. This regime is indicated as 'Complex conjugate roots' in figure 6, and can only occur if both Bi 3.87 and Pe > 18.9. Most physical studies involve small values of Bi, and so this condition is not particularly restrictive. We note that none of the other models have any comparable positivity results; the integrands would involve linear combinations of Θ and its derivatives with respect to x, and could not be bounded using only conditions on the sign of Θ.
The derivation of our new model supplies an approximation of the interior temperature profile in addition to an equation for the surface temperature. Contour plots of this interior temperature are shown in figure 8 for various Pe with Bi = 0 in response to strip heating. The predicted interior temperature profiles appear broadly similar to the full 2-D results. The biggest deviations are for Pe = 4 and Pe = 6, which are closest to the critical value of Pe * = 4.592. Even for these values of Pe, the contours are in very good agreement away from the regions of rapid change.
Our final thermodynamic consistency check comes from property P4; there should be no extrema of the temperature inside the fluid layer. For the cases shown in figure 8, this property is satisfied. However, if we repeat the calculations from .592, so this case is close to the change in criticality. Although S(x) 0 here, the internal temperature field has a region where θ < 0 (inside the blue contours near x = −2), and also a local maximum (inside the red contour near x = 2).
figure 8 with Pe close to the critical value Pe * , we can construct cases with negative internal temperature. In figure 9, we show a contour plot of the internal temperature for one such case, with Bi = 0 and Pe = 4.5. Although the surface temperature is non-negative, there is a region of negative temperature inside the layer just upstream of the heating strip and also a local maximum of θ inside the fluid layer above the downstream end of the heating strip.

Time-periodic heating
We have shown that our new model (4.19) successfully predicts steady heat transfer through a uniform flow, at least for cases that are not too near to the critical Péclet number. However, time-dependent heating is also important, particularly for applications involving real-time control. The simplest test of time dependence is to take the applied heating and resulting temperature profile to be sinusoidal in time. For the long-wave models, we can obtain solutions straightforwardly using Fourier transforms in x. Fourier transforms are also appropriate for the full 2-D system, though we must numerically solve an ordinary differential equation (ODE) in y at each k. The results of this calculation for a periodically oscillating strip are shown in figure 10, and show very good agreement between the projected and full system for these parameter values. Note that, in a full initial-value problem, the periodic response should be combined with a full start-up problem, but the decay of transients in this initial problem has been addressed in § 5.1.

Discussion
With respect to the five thermodynamic properties introduced in § 3, the new projected model satisfies P1-P3, and satisfies P5 (unless Bi > 3.87), but can fail P4 for some temperature profiles if Pe is close to the critical value. The performance of all the long-wave models against the five properties is summarized in table 1. The explicit (Benney) model fails catastrophically once used outside the long-wave limit. The  6. For the other models, at any Pe and Bi we can construct a temperature profile that would yield a negative surface temperature.
weighted-residual and parabolic models both make much more sensible predictions, but, because they rely on derivatives of the wall temperature to obtain accurate results, they are not robust when applied to heterogeneous heating with short-wave components. Our new projected model is much better behaved, and indeed always obeys P1-P3, but suffers from a singularity resulting from the change in direction of one of the characteristics from upstream to downstream as the Péclet number increases. This means that, near the critical value of Pe, we can construct examples that produce negative internal temperature. Hence our conclusion from the analysis of second-order long-wave models (in the case of Ma = 0) is that no such model can always obey properties P1-P3 and P5 for arbitrary substrate heating while making meaningful predictions about the interior temperature profile that obey property P4. Nonetheless, the new projected model is more robust than the alternatives and hence is a valuable new tool to analyse the behaviour of falling liquid films in response to variable heating.

Coupling to film equations
The fluid temperature affects the hydrodynamic problem through tangential stresses at the air-fluid surface, and hence only via S(x, t); this coupling to hydrodynamics has been well studied within the context of systems with uniformly heated walls. On the other hand, the heat transfer problem is influenced by the velocity field throughout the fluid layer, and this coupling can become quite intricate. Our aim here is to derive long-wave momentum and heat equations both valid for non-uniform heating.

Effect of surface temperature on film equations
We use a two-equation, weighted-residual model for the film dynamics, coupling the film height h(x, t) and the depth-integrated flux q(x, t), which naturally satisfy where the second equation follows from mass conservation. We follow the standard weighted-residual approach, seeking to include surface tension, inertia, axial diffusion and Marangoni stresses. In order to retain second-order terms relating to Marangoni effects, and including some interactions between inertia and Marangoni stresses, we set Re = O(δ), Ma = O(δ −1 ), Ca = O(δ 2 ) and reach In this case the Marangoni number appears both through tangential stresses and as variable surface tension from the normal stress balance.

Heat transfer through a non-uniform film
The heat equation is affected by film non-uniformity in two ways: the boundary conditions should be applied at y = h(x, t), and the velocity field (u, v) appears directly in the heat equation. The first effect here is easy to deal with, but, for the second, we need a suitable characterization of the velocity field. Noting that the derivation of the equation for surface temperature makes use of u and its derivatives at the surface, it is reasonable to use the flux-parametrized velocity profile, as is also used in the derivation of the weighted-residual hydrodynamic equations: This velocity field is incompressible, and obeys u y = 0 and h T = −q X =v − uh X at y = h. As for the uniform film case, we derive the projected heat equation by posing an expansion of θ (X, y, T) for small δ in the form (4.4). The boundary conditions on θ 0 , θ 1 and θ 2 are obtained from expanding the two conditions θ = S(X, T), θ y − δ 2 h X θ X = −Bi S(1 + δ 2 h 2 X ) 1/2 , on y = h(X, T), (6.4a,b) in powers of δ. Similarly, the PDEs obeyed by θ 0 , θ 1 and θ 2 are determined from a small-δ expansion of δPe(θ T + uθ X +vθ y ) = δ 2 θ XX + θ YY . (6.5) We close the model with the condition that θ (X, 0, T) = Θ(X, T). The resulting equation involves time derivatives of the hydrodynamic variables h and q. We can use h T = −q X to eliminate time derivatives of h, to reach (in terms of the original variables) We can also include the effect of Marangoni stresses on the flow field used to derive the heat equation: For steady states, and assuming that Ma and Pe are both O(1), the Marangoni stresses have the effect of adding the term (6.9) to the left-hand side of (6.7). (The expression required to describe Marangoni corrections to the unsteady equation (6.6) is considerably longer.) 7. Steady states subject to localized heating 7.1. Small displacement solutions In the limit of small Ma, we can obtain a solution by linearizing about a flat film. We let We can analyse the linear solution via a heat transfer problem, which relates Θ and S, and a hydrodynamic problem, which relates S and h. We solve the hydrodynamic problem (6.2) at O(Ma) to obtain The heat transfer problem, to determineŜ, is exactly the linear problem discussed in § 5. For small displacements of the fluid surface, we can use an Orr-Sommerfeld analysis to solve both the Navier-Stokes and heat equations and hence determine the equivalent full 2-D surface displacement. The results are shown in figure 11 for comparison with the long-wave results, in a case with strong surface tension and low Reynolds number, so that the base flow is stable. The long-wave results are in very good agreement with those from the linearized 2-D system. Note that the half-width of the strip heating profile has been increased from X = 2 to X = 5; this leads to increased deflection of h(x) and makes the surface shape easier to visualize.

Finite-amplitude solutions
At small Ma, the heat problem is exactly the same as that for flat films. Hence, to test the validity of the dependence of (6.7) on h and q, we must increase the deformation, most easily achieved by increasing Ma. For large-amplitude solutions, we solve numerically in a periodic domain, applying the strip-like distribution (5.2) for Θ and obtain steady states by continuation in Ma from h = 1 at Ma = 0, subject to the constraint that the mean layer height is fixed at 1. Results using the long-wave equations (6.2) and (6.7), computed using the bifurcation and continuation transfer package AUTO 07P (Doedel & Oldman 2009), are shown in figure 12. With all other parameters fixed, we find that increasing Ma increases the surface deformation h, accompanied by only small changes in the surface temperature profile S. Solutions are available for a wide range of Ma (at least up to Ma = 25), except for the two cases with Pe = 5. This value of Pe = 5 is close to the critical Pe for a uniform film with h = 1, q = 2/3 for both Bi = 0 (where Pe * = 4.592) and Bi = 1 (where Pe * = 4.917), so that the coefficient of S xx , at h = 1 and q = 2/3, is small but non-zero in both cases. Therefore, it is possible that the coefficient of S xx in the nonlinear calculations, 3) may vanish as q and h vary. For Bi = 0, the sign of (7.3) is a function only of Pe and q, and would vanish at q = 0.6123 for Pe = 5, for example. As q is spatially constant in steady calculations, 874 A. B. Thompson, S. N. Gomes, F. Denner, M. C. Dallaston and S. Kalliadasis in fact this would mean that the coefficient of S xx vanishes everywhere in the domain simultaneously if q exactly coincides with the critical value, but the equation would be non-singular for all other q. In the case Pe = 5, Bi = 0, continuation from Ma = 0 yields solutions up to Ma = 14.68, where q = 0.614208. In order to obtain solutions for larger Ma, we use continuation from Pe = 1, Bi = 0, first increasing Pe and then decreasing Ma. Using this method, the AUTO code yields converged solutions down to Ma = 14.88.
In contrast, for Bi = 0, the dependence of the sign of the coefficient of (7.3) on h means that the coefficient can vanish at some points in the domain but not others, and hence, for sufficiently large surface deformation, we expect the ODE (6.7) to have one or more singular points. Fixing Pe = 5, Bi = 1 and q = 2/3, (7.3) would vanish at points where h = 1.326. However, the calculations shown in figure 12 are performed with a constraint on layer height, so that q decreases slightly as Ma increases. For Bi = 1, Pe = 5, the last converged solution when increasing Ma from zero is for Ma = 7.7899, where q = 0.66237 and the maximum value of h is 1.15829. Using a continuation method from solutions at large Ma for Pe = 1, we can also obtain solutions at Pe = 5, Bi = 1 for Ma > 17.64, where q = 0.64321 and the minimum value of h is 0.684856. Comparison of these values of Pe, Bi, h and q with (7.3) suggests that the border of the region where we can compute solutions does indeed correspond to conditions where the coefficient of S xx locally vanishes.
The finite-amplitude calculations require 2-D Navier-Stokes simulations for validation. Steady solutions are calculated using the finite-element library oomph-lib (Heil & Hazel 2006). In this computation, Marangoni effects are created by specifying the variation of surface tension with temperature, leading to effects in both the normal and tangential components of the dynamic boundary condition, consistent with the boundary conditions (2.16) and (2.17) used to derive our long-wave equations. The air layer is treated as a region of constant pressure and temperature, and heat loss is parametrized by the same Biot-number condition as we used to derive the model. The Navier-Stokes calculations are performed in closed conditions: a periodic domain with mean layer height 1.
Several of the long-wave calculations are in remarkably good quantitative agreement with the full nonlinear Navier-Stokes results. In particular, for Bi = 1, the long-wave results shown in figure 12 are in excellent agreement for h(x) across the range of Ma where they exist. The predictions for S(x) are also in good agreement at Bi = 1, capturing the notable increase in S at the downstream end of the heater as Ma increases. The agreement for Bi = 0 is less good than for Bi = 1; this is likely due to the generally reduced surface temperatures for larger Bi leading to smaller variations in h(x) at the same Ma and hence less significant nonlinearity. For Ma 10 and Bi = 0 at the larger Pe, the predictions for the upstream fluid ridge have the wrong shape when compared to the 2-D nonlinear results. This may be related to the observation that, while the long-wave results correctly predict an increasingly rapid decay of S(x) downstream of the heater as Ma increases, they do not capture the slowing of the upstream decay, which should lead to increasingly important upstream influence regions.
Similar calculations including the Marangoni influences on the heat equation itself via (6.9) are shown in figure 13. For Pe = 1, Bi = 0, these correction terms lead to better agreement for S(x), particularly just upstream of the heating strip. There is also slightly improved agreement for h(x). However, for larger Pe, the inclusion of Pe Ma terms via (6.9) reduces the range of Ma over which solutions can be obtained. In particular, note that for Bi = 0, we obtain solutions at Pe = 5 up to Ma = 9.39 and at Pe = 10 up to Ma = 18.0. In both cases, the last converged solution is associated with strong deviations from the Navier-Stokes results for S(x) at the upstream edge of the heating strip. The full 2-D solutions persist without qualitative changes in behaviour until at least Ma = 25. We can isolate the effects of Pe Ma interactions in a much simpler problem: the zero-Ca, zero-Re limit of strong surface tension, where the normal component of the dynamic boundary condition implies simply that the surface is flat, but Marangoni stresses drive an internal flow due to the tangential stress balance. If S(x) is known, the full equations for u and v are linear: (7.4a−c) subject to u = v = 0 on y = 0, and u y = −Ma S x , v = 0 on y = 1. The unique solution for u and v with S given can be obtained by Fourier transforms. However, the coupled problem is nonlinear, as the surface temperature S(x) = θ (x, 1) is determined by the heat equation, which itself involves the velocity field Within the long-wave limit, the structure of the equations for this limiting case is similar to that of the full problem. The steady velocity field induced by a surface temperature S(x) is given by It is straightforward to verify that the 'projected' heat equation can also be obtained by setting h = 1 and q = 2/3 in the equivalent non-uniform equation (6.7) with the addition of (6.9): Pe 2 239 10 080 + Without the Pe Ma term, the velocity field is independent of S, and so if h is known, u and v are also known, the heat equation is linear, and the full and long-wave heat equations each have a unique solution for S. However, with non-zero Pe Ma, there is no reason to expect a unique solution in either the full or long-wave system. We can apply the same strip heating profile as in figure 13 and calculate the surface temperature profile in the strong surface tension limit. The long-wave and full 2-D results (calculated using the full Navier-Stokes code with Ca = 10 −5 ) are shown in figure 14. We find that the long-wave solution branch for Pe = 10, Bi = 0 persists to Ma ≈ 40, but the maximum value of S increases with increasing Ma and in fact crosses 1 at Ma ≈ 32, thus violating property P5, which should still hold in this limit. For relatively small Ma (up to Ma = 20 for Bi = 0 and up to Ma = 40 for Bi = 1), the full 2-D results are quite similar to the long-wave results, with a broadening region of constant temperature above the strip. Beyond these Marangoni numbers, the 2-D results predict that this region continues to broaden, and the surface is always hottest above the downstream end of the strip. At the very largest Ma (Ma = 100 for Bi = 1), we start to see the appearance of a reversal in the direction of the temperature gradient, in qualitative agreement with the long-wave results. However, the long-wave model predicts that this feature should emerge at much lower Ma and be much larger in magnitude by Ma = 100. The long-wave models hence lose accuracy at around Ma = 20, with Pe = 10, by which point the combination Pe Ma is clearly not small. Hence, if better agreement is required at these parameters, the equations should be rederived under the assumption of large Pe Ma with respect to the long-wave parameter δ, likely introducing further nonlinear interactions.

Conclusion
We have explored in detail the formulation of long-wave models for predicting the effect of localized wall heating on surface temperature, and via variation of surface tension with temperature, on the dynamics of falling liquid films. This problem can be split into two components: the relationship between surface temperature S and wall temperature Θ, which is moderated by flow, and the relationship between surface temperature and film dynamics. The latter coupling occurs via Marangoni stresses in the momentum balance, and has already been studied in a variety of contexts. However, the first part, on the transmission of localized heating through a uniform or non-uniform film, is not yet well understood and our paper represents the first comprehensive study in this area.
We began our analysis of the heat transmission problem by deriving thermodynamic constraints for the full 2-D heat equation subject to a uniform Nusselt flow. We proved five properties: (P1) the magnitude of temperature variations at the fluid surface must be less than that imposed at the wall; (P2) at short wavelength, the surface feels no temperature variation; (P3) in the absence of imposed heating, the film should decay to uniform temperature; (P4) there can be no local extrema of the temperature field inside the fluid layer; and (P5) the surface temperature must be bounded by the extrema of the wall and air temperature.
In previous studies of film patterning by imposed steady heating, most authors have chosen to use only the leading-order long-wave model S = Θ/(1 + Bi h). This leading-order model takes no account of advection or streamwise diffusion, and predicts instantaneous reaction to any change in wall temperature. The leading-order model is especially inaccurate in the large-Pe regime that would apply if the film Reynolds number is large enough for the flow to be unstable.
In many isothermal falling liquid film studies, the hydrodynamic model used is second-order-accurate or higher with respect to the long-wave parameter, and so we focused here on thermal models that are also second-order-accurate. If the film is uniform (always attainable with Ma = 0), the surface temperature equation should be linear in S and Θ, and there are five free parameters in choosing a second-order-accurate model with at most second-order derivatives of S and Θ with respect to space and time. The most obvious of these models is a Benney-type model, yielding the surface temperature explicitly in terms of Θ, Θ x , Θ t , Θ xx , Θ xt and Θ tt . However, although this model is quantitatively accurate in the long-wave limit, it is divergent when evaluated in the short-wave limit, gives completely inaccurate results when applied to a smoothed heating strip, and fails all applicable thermodynamic tests. For the purely hydrodynamic isothermal problem, the deficiencies of the Benney model can be remedied by use of a weighted-residual model, which leads to an additional degree of freedom in the system. However, application of the same methodology for the heat equation is less successful; although it indeed leads to an evolution equation for S rather than an explicit expression, it still fails most of the thermodynamic tests because it still relies upon derivatives of Θ, which is not even required to be smooth in the full 2-D heat equation.
We are therefore led to seek a second-order model that gives S implicitly, but with only Θ and not its derivatives appearing as a source function. There is a unique model of this form among the class of linear second-order-accurate models. We showed that this model can be obtained systematically by solving for a temperature field that exactly satisfies the two boundary conditions at the surface along with the bulk equation. We can then 'read off' the wall temperature in terms of S and its derivatives. Requiring that this wall temperature coincides with the imposed function leads to an evolution equation for S. The resulting model always satisfies P1-P3.
Having focused our attention on this new 'projected' model, we proceeded to examine its behaviour across a wider range of parameters, for both uniform and non-uniform films. For uniform films, the equation is always hyperbolic, but the direction of its characteristics depends on Pe. For Pe below a critical value Pe * , one family of characteristics travels upstream and the other downstream and, for both steady and time-dependent systems, a localized heating strip affects the surface temperature both upstream and downstream. For Pe above the critical value, both families of characteristics travel downstream, and a heating strip has strictly downstream influence. The value of Pe * depends slightly on Bi, ranging from 4.592 to 6.563 as Bi ranges from zero to infinity. If the film is uniform and the system is steady, the projected model reduces to a second-order, constant-coefficient ODE for the steady surface temperature corresponding to a given wall profile. Furthermore, solving this ODE via variation of parameters allows us to show that, except in a subset of the region Bi > 3.87, non-negative wall temperature leads to non-negative surface temperature; thus P5 holds.
Our new model is unique among second-order models in satisfying P1-P3 and P5, but still does not guarantee that the corresponding internal temperature profile (supplied as part of the derivation) is thermodynamically consistent with property P4. By considering cases with sharp wall temperature gradients, with Pe close to Pe * , we can construct counterexamples to P4. The question of negative surface and internal temperatures also arises in the study of travelling waves on films subject to uniform heating, and has not been conclusively avoided there either. However, we should also note that the constraints on modelling imposed by heterogeneous heating do not apply for uniform heating, so that a wider class of models could be relevant in that case.
Having examined the performance of our new model against the five thermodynamic properties, we next considered how film non-uniformity modifies the temperature equation. The derivation follows in a similar way to that for a uniform film, but it requires explicit parametrization of the velocity field, and can in general lead to a rather complicated governing equation. If h and q are predetermined, the evolution equation is again hyperbolic, but, as the coefficients in the long-wave heat equation depend on h and q, the direction of the characteristics may evolve with the solution.
As the Marangoni number is increased for the case of steady heating, our model predicts increasing surface deflection, in accordance with full Navier-Stokes results. The variations in h(x) are accompanied by minor changes in S(x), at least up to Ma ≈ 20. There is no evidence of bifurcations at smaller Ma, though calculations for Pe near to Pe * do not converge for all Ma, due to the development of local singularities in the equations associated with the change in criticality of the equations.
At larger Ma, the full 2-D results still show only mild changes in S(x). Some of these features are well predicted by the long-wave models; for example, the increase in surface temperature at the downstream end of the heater for non-zero Bi. We can seek to improve the agreement by including Marangoni corrections in the velocity field used to derive the heat equation; a nonlinear effect associated with the combination Pe Ma. At small Pe, these corrections indeed lead to better agreement, but for Pe = 10 and Ma ≈ 20, the long-wave model yields unphysical predictions such as an increasingly hot free-surface near the upstream end of the heater. This upstream heating effect emerges even if film non-uniformity is eliminated by considering the limit of very large surface tension, where the film surface should be flat, and can lead to very hot regions on the air-fluid surface, violating P5. These predictions can likely be improved to some extent by allowing Pe Ma to be asymptotically large in the derivation of the long-wave model, but such broad applicability would come at the cost of further increase to the model complexity.
In our analysis we have assumed throughout that the temperature at the wall is exactly specified. However, in physical implementation it is more likely that the boundary condition would involve either a specification of the heat flux or some linear combination of temperature and heat flux (conceivably coupled to another thermal problem within the wall such as considered by Dallaston, Tseluiko & Kalliadasis (2016)). Our modelling methodology can be easily extended to these more general boundary conditions (see appendix A). The derivation again involves expanding the temperature field about the surface values; the model is closed by the condition that the surface-based expansion for θ be consistent with whatever boundary conditions are required at the wall. We note that the analysis would still hold even if parameters in the boundary condition, for example, the thermal conductivity of the wall, are themselves variable. This would also allow long-wave access to boundary conditions of mixed type, and more accurately reflect the experimental set-up for the analysis of ridge instability, where a heating strip of known heat flux is embedded in either an insulating wall, or one kept at a constant temperature away from the heater.
The long-wave modelling methodology developed here should serve as a useful modelling framework for a number of problems, including our motivating problem of control via spatially and temporally varying heating, and also the canonical system of the spanwise instabilities of the capillary ridge induced by a single heating strip. Our analysis can also be modified to account for a nonlinear dependence of surface tension on temperature. However, we have not considered the case of temperature-dependent viscosity, which remains a significant challenge.