Solitary waves on falling liquid films in the inertia-dominated regime

We offer new insights and results on the hydrodynamics of solitary waves on inertia-dominated falling liquid films using a combination of experimental measurements, direct numerical simulations (DNS) and low-dimensional (LD) modelling. The DNS are shown to be in very good agreement with experimental measurements in terms of the main wave characteristics and velocity profiles over the entire range of investigated Reynolds numbers. And, surprisingly, the LD model is found to predict accurately the film height even for inertia-dominated films with high Reynolds numbers. Based on a detailed analysis of the flow field within the liquid film, the hydrodynamic mechanism responsible for a constant, or even reducing, maximum film height when the Reynolds number increases above a critical value is identified, and reasons why no flow reversal is observed underneath the wave trough above a critical Reynolds number are proposed. The saturation of the maximum film height is shown to be linked to a reduced effective inertia acting on the solitary waves as a result of flow recirculation in the main wave hump and in the moving frame of reference. Nevertheless, the velocity profile at the crest of the solitary waves remains parabolic and self-similar even after the onset of flow recirculation. The upper limit of the Reynolds number with respect to flow reversal is primarily the result of steeper solitary waves at high Reynolds numbers, which leads to larger streamwise pressure gradients that counter flow reversal. Our results should be of interest in the optimisation of the heat and mass transport characteristics of falling liquid films and can also serve as a benchmark for future model development.


492
F. Denner and others

Introduction
Falling liquid films have been a topic of fundamental and applied research for several decades since the pioneering experiments by Kapitza (1948) and Kapitza & Kapitza (1949). A falling liquid film is a convectively unstable open-flow hydrodynamic system with a rich variety of spatiotemporal structures and a sequence of wave instabilities and transitions that are generic to a large class of hydrodynamic and other nonlinear systems. Despite their apparent complexity one can still identify robust coherent structures, i.e. solitary waves, in what appears to be a randomly disturbed surface. These structures interact continuously with each other as quasiparticles. Hence, a falling liquid film can serve as a canonical reference system for the study of spatiotemporal chaos and weak/dissipative turbulence. In addition, falling liquid films are typically associated with low flow rates, low pressure drops, small thermal resistances and large contact areas per unit volume. Not surprisingly, they play a central role in a wide spectrum of engineering applications. Prominent examples of applications involving falling liquid films are evaporators-condensers, heat exchangers, chemical reactor columns and cooling devices. Falling liquid films are also exploited in saltwater desalination plants, in the food processing industry as well as in coating processes. Furthermore, falling liquid films are found in problems related to earth science and geophysics, for instance underwater gravity currents, lava flows and the propagation of bores in rivers. Therefore, an accurate understanding of the hydrodynamics of falling liquid films is crucial in the design and optimisation of processes and devices that exploit these flows. The present study focuses on the hydrodynamic mechanisms that govern solitary waves on inertia-dominated falling films and, in particular, on flow features that have been shown to improve the heat and mass transport characteristics of the films.

Theory and physical aspects
A falling liquid film is unstable to long-wave perturbations above a critical Reynolds number (Yih 1963) due to an instability mechanism that is mainly driven by shear stresses acting at the interface (Hinch 1984;Kelly et al. 1989;Smith 1990). In general, inertia and the streamwise component of gravity destabilise the flow, whereas surface tension and the cross-stream component of gravity have a stabilising effect. A natural source of perturbations to a falling film flow is white noise at the inlet, whereby the resulting instability selects the fastest growing mode of the system followed by a secondary modulation instability that converts the primary wave field into solitary waves (Liu, Paul & Gollub 1993;Chang 1994;Kalliadasis et al. 2012). Liu & Gollub (1994), however, reported regular waves by imposing a monochromatic perturbation at the inlet in their experimental arrangement. In this case, the instability synchronises with the inlet forcing frequency Chang 1994), provided that the forcing amplitude is sufficiently large.
At low forcing frequencies, the long-wave instability develops into fast-moving solitary waves that consist of a primary wave hump preceded by capillary ripples at the front. These waves are predominantly two-dimensional (2D) structures (Leontidis et al. 2010), and in the inertia-dominated regime (typically for Reynolds numbers, Re 20-200) the velocity profile can deviate strongly from its originally parabolic shape (Malamataris, Vlachogiannis & Bontozoglou 2002;Scheid, Ruyer-Quil & Manneville 2006). In addition, if the inertia of the liquid film is sufficiently high, the maximum velocity of the fluid exceeds the wave speed of the solitary waves, leading to a recirculation zone in the moving frame (Maron, Brauner & Hewitt 1989).
Solitary waves on falling liquid films in the inertia-dominated regime 493 This zone, which appears in the main hump of the solitary waves, can considerably impact the heat and mass transport in the film due to the enhanced mixing (Roberts & Chang 2000;Albert, Marschall & Bothe 2014). Previous studies have found that the wave frequency, the inclination angle of the substrate, as well as viscous dispersion and capillary effects have a significant influence on the extent of the recirculation zone (Miyara 2000;Rohlfs & Scheid 2015).
The number of capillary ripples preceding a solitary hump increases with decreasing forcing frequency, increasing Reynolds number and decreasing viscous dispersion (Nosoko et al. 1996;Pradas, Tseluiko & Kalliadasis 2011), in order for surface tension to balance inertia. Capillary ripples propagate at the same speed as the main hump of the solitary wave (Dietze 2016), since they are continuously compressed by the main hump. As a result, faster capillary ripples have a larger amplitude (Dietze 2016), in contrast to freely oscillating capillary waves which exhibit a smaller phase velocity for increasing initial wave amplitude (Denner, Pare & Zaleski 2017a). Due to the high interface curvature of the capillary ripples, pressure gradients in the streamwise direction are significant and may lead to flow separation and local reversal of the flow direction in the laboratory frame of reference. This phenomenon was first noted by Portalski (1964) and Massot, Irani & Lightfoot (1966), and was later investigated in more detail in terms of the underlying mechanisms (Dietze, Leefken & Kneer 2008;Dietze, Al-Sibai & Kneer 2009) and the influence of inertia, wave frequency and inclination angle (Malamataris & Balakotaiah 2008;Chakraborty et al. 2014;Rohlfs & Scheid 2015). Numerical simulations conducted by Kunugi & Kino (2005) demonstrated a considerable enhancement of the heat and mass transport characteristics of the film as a result of flow reversal.

Experimental and numerical methods
From an experimental point of view, a range of optical (typically laser-based) methods have been developed for the investigation of the instantaneous motion and spatial extent of unsteady interfacial flows and, specifically, liquid film flows (Lel et al. 2005;Schlagen et al. 2006). Fluorescence intensity-based measurements of the instantaneous and local film height have been conducted, for example, in laminar falling films , in annular gas-liquid flows (Alekseenko et al. 2009(Alekseenko et al. , 2012 and in gas-sheared liquid flows in a horizontal rectangular duct (Cherdantsev, Hann & Azzopardi 2014). Liu and co-workers ) determined the critical Re for the onset of wave formation as a function of the inclination angle, and a phase boundary which separates the saturated, finite-amplitude wave regime from the multipeaked solitary wave regime was identified. In the work of Alekseenko et al. (2009Alekseenko et al. ( , 2012 and Cherdantsev et al. (2014), 2D and three-dimensional (3D) film-thickness measurements were used to scrutinise the topology of large-amplitude disturbance waves as a function of Re, and to link the generation of ripples forming atop these waves to the mechanisms of liquid entrainment. Markides and co-workers (Zadrazil, Matar & Markides 2014a; used planar laser-induced fluorescence imaging (PLIF) instead, in order to investigate the coupling between the film-thickness and velocity-field variations over a 2D domain. In Mathie, Nakamura & Markides (2013) and Markides, Mathie & Charogiannis (2015), the same technique was employed alongside infrared thermography (IR), thus allowing for interface-temperature measurements to be retrieved simultaneously with film-thickness measurements. Based on this approach, the authors reported significant heat transport enhancements depending on film topology.

F. Denner and others
Optical velocimetry techniques can broadly be categorised depending on whether a molecular-tagging or particle-seeding approach is adopted. The photochromic dye activation (PDA) technique stands as a prime example of the former, allowing for direct velocity-profile measurements via the generation of thin fluorescing traces arranged perpendicular to the direction of the flow (Karimi & Kawaji 2000;Moran, Inumaru & Kawaji 2002). Particle-image velocimetry (PIV) and particle tracking velocimetry (PTV), instead, allow for 2D as well as 3D velocity-field characterisation by tracking the motion of seeded particle groups (PIV) or individual particles (PTV). Examples of the use of micro-PIV include the study of laminar, vertically falling films (Adomeit & Renz 2000) and the capillary-wave region (Dietze et al. 2008), while PIV/PTV has been used, for example, in the study of downward co-current gas-liquid annular flows (Zadrazil, Matar & Markides 2014b), where the presence of multiple recirculation zones within the disturbance waves was revealed.
In terms of modelling and computational methodologies, a substantial number of studies have focused on low-dimensional (LD) formulations, which simplify the governing equations and associated wall and free-surface boundary conditions considerably. Using the long-wave expansion, a hierarchy of model equations has been obtained, starting from the so-called Benney equation (Benney 1966) whose region of validity is restricted to near-critical conditions (Kalliadasis et al. 2012), to integral-boundary-layer equations which are valid away from criticality and describe the evolution of both the film height h and the local flow rate q (Shkadov 1967) (see also the work by Scheid et al. (2006), where different LD models are compared to each other). Numerical modelling has also received considerable attention, with direct numerical simulation (DNS) finding increasing application in the study of falling liquid films. The work of Ramaswamy, Chippada & Joo (1996) was the first DNS study to conduct a detailed investigation of this problem, focusing on the stability of falling films as well as on the development and interaction of interfacial instabilities. Subsequently, Malamataris et al. (2002) showed that the velocity profile of the liquid film is strongly non-parabolic in solitary waves, whereas studies of Gao, Morley & Dhir (2003) and Nosoko & Miyara (2004) investigated the dynamics of periodically perturbed falling films at different forcing frequencies and Reynolds numbers. A common assumption made in the majority of computational studies is a periodic computational domain in which the film thickness is conserved, hence increasing the flow rate at the onset of interface waves (Ramaswamy et al. 1996;Gao et al. 2003) and imposing a particular wavelength on the instabilities (Kalliadasis et al. 2012). In reality, however, the system is not closed but open and, therefore, it is not the mean film thickness but the mean flow rate that is conserved, resulting in a decrease of the average film thickness for wavy film flows (Chang, Demekhin & Kalaidin 1995;Gao et al. 2003). Numerical studies using periodic domains are, hence, limited in their ability to accurately predict spatiotemporally evolving interface instabilities.  (Denner et al. 2016), a yet unidentified mechanism must be responsible for the saturation (stagnation or reduction) of the height of solitary waves. Similarly, with respect to the flow reversal observed underneath the trough preceding the solitary wave, the hydrodynamic mechanisms responsible for the upper limit in Reynolds number above which no flow reversal was observed in simulations of Chakraborty et al. (2014) have yet to be explained. It is also worth noting that none of the previous studies included a detailed and systematic comparison between experimental measurements and numerical results using LD models or DNS. If such a comparison confirms good agreement, a concurrent application of these methods may provide more detailed insight into the complex hydrodynamics of solitary waves.
Aiming to provide answers to these open questions, we adopt a synergistic approach with a combination of experiments, DNS as well as LD modelling. A direct comparison of experimental and numerical results is instrumental in assessing the predictive quality of the available computational methods and in analysing the flow field and wave characteristics in detail. We find that experimental measurements and DNS results agree very well, in particular with respect to film height and flow velocity at the interface, while the applied LD model is favourably compared with both of them. Of particular interest in our study is the flow field within the falling liquid films and its influence on the dynamics of solitary waves and capillary ripples, in order to elucidate the mechanisms stabilising the wave height and influencing the flow separation and reversal at large Reynolds numbers. The presented experimental and numerical results can also serve as a benchmark for future studies of the hydrodynamics of falling liquid films.
In § 2 the characterisation of falling liquid films and the associated solitary waves is discussed. In § 3 we introduce the experimental and numerical methods we adopt and detail the case studies we undertake. The results are presented and discussed in § 4 and conclusions are offered in § 5.

Characterisation of falling liquid films
Consider a film flowing down a planar inclined plane forming an angle β with the horizontal direction, as shown in figure 1. There exists a spatially uniform solution, often referred to as the Nusselt flat-film solution (Nusselt 1916), with film height h N = 3 3µ l q N ρ l g sin β (2.1) and average film velocity where µ l and ρ l are the dynamic viscosity and the density of the liquid film, respectively, q N = u N h N is the flow rate per unit span and g is the gravitational acceleration. Based on this idealised solution, various characteristic scalings can be formulated.
2.1. Nusselt scaling Making use of the representative pressure scales of the three major hydrodynamic mechanisms acting on the falling film (inertia, viscous stresses and surface tension) h(x, t) g FIGURE 1. Sketch of the profile geometry for a liquid film on a substrate with inclination angle β to the horizontal. h(x, t) is the local film thickness with respect to a Cartesian coordinate system (x, y), with x the streamwise coordinate and y the outward-pointing coordinate normal to the substrate.
several dimensionless parameters can be defined. The Reynolds number expresses the relative importance of inertia over viscous effects, the Weber number where σ is the surface tension, represents the relative importance of inertia to surface tension, and the capillary number, a measure of the relative importance of viscous stresses and surface tension. The Reynolds number Re can also be regarded as the dimensionless flow rate. These non-dimensional numbers are related to each other by We = Re Ca and are often referred to as Nusselt scaling in the context of falling liquid films, because they are based directly on the Nusselt flat-film solution. Numerical results presented by Denner et al. (2016) suggest that the capillary number is not suitable to characterise inertia-dominated solitary waves. Denner et al. (2016) recently proposed an 'inertia-corrected scaling' for the characterisation of solitary waves on falling liquid films based on the driving Nusselt velocity u * N = u N sin β. Using this velocity scale, the Reynolds and Weber numbers are redefined as Re * = Re sin β and We * = We(sin β) 2 , respectively, which was shown to lead to a self-similar characterisation of the shape and dispersion of solitary waves on falling water films for a given inclination angle of the substrate and forcing frequency.
Unlike the dimensionless numbers above, which depend on the flow rate, the Kapitza number where ν l = µ l /ρ l is the liquid kinematic viscosity, depends only on the properties of the liquid and inclination angle β. Hence, for a fixed liquid and β, the Kapitza number is fixed and the only free parameter is the Reynolds number. can be obtained. This reduced Reynolds number δ compares inertia and viscous effects at the length scale imposed by surface tension and gravity. On the other hand, the viscous dispersion number η appears along with every second-order streamwise viscous term in the momentum equation and has a dispersive effect on the speed of linear waves (Pradas et al. 2011;Kalliadasis et al. 2012). The Shkadov scaling is particularly useful for the characterisation of LD models, as further discussed with respect to the applied LD model in § 3.3. Based on the reduced Reynolds number δ, Ooshida (1999) identified two regimes: the drag-gravity regime for δ 1 and the drag-inertia regime for δ 1. In the draggravity regime, the flow is similar to that of the Nusselt flat-film solution, with both inertia and surface tension playing effectively a perturbative role (Kalliadasis et al. 2012). In the drag-inertia regime on the other hand, inertia and surface tension are no longer corrections to the Nusselt flow, the cross-stream velocity component becomes significant and the velocity profile of the film flow can become strongly non-parabolic (Malamataris et al. 2002). Recent studies (Chakraborty et al. 2014;Rohlfs & Scheid 2015) found the transition between the drag-gravity and drag-inertia regimes to occur at 1 δ 6. Additional effects -for example, shear-thinning and non-Newtonian behaviour -can lead to complex behaviour in the vicinity of this transition (such as hysteresis) (Pradas et al. 2014).

Experimental methodology
The experimental apparatus and methodologies developed to generate the film thickness and velocity data have been described extensively in our previous publications (Charogiannis et al. , 2017, and therefore only a summary is presented here for the sake of clarity and completeness. A schematic of the experimental arrangement, itself an evolution of a previously developed set-up (Mathie et al. 2013;Markides et al. 2015), is depicted in figure 2. The apparatus consists of a main test section over which film flows develop, and a closed loop via which the liquid circulates. Prior to entering a flow-preparation distribution box, the flow is split into a steady supply Q and a pulsating supply Q . The latter is generated by an oscillator valve that allows accurate control of the forcing frequency and amplitude to be achieved, selectively triggering the rapid growth of fully developed wave regimes within the confines of our test section. The instantaneous flow rate into the distribution box is calculated from the output signal of an ultrasonic flowmeter. The role of the distribution box is to create a well-characterised uniform flow with low noise levels, and to introduce this over the test section with uniform thickness.

Test section and imaging arrangement
The test section comprises a thin soda-lime glass plate with 0.7 mm thickness, mounted on an aluminium support frame and inclined at 20 • to the horizontal. The Reynolds number can be evaluated at the flow inlet using the mean flow-rate measurement per unit width of the channel Q, as where U box is the bulk velocity inside the box and D is the channel depth. A frequency-doubled Nd:YAG laser (532 nm) is used to illuminate the flow, which is seeded with Rhodamine B dye for PLIF and glass-hollow-sphere particles for PIV, at a sampling rate of 100 s −1 . Sheet optics are employed in forming a thin (approximately 0.1 mm) laser sheet extending along the (streamwise) length of the substrate, while a pair of LaVision HS 500 CMOS cameras equipped with Sigma 105 mm f/2.8 Macro lenses are employed to achieve the desired magnification (between 28.0 µm and 29.7 µm per pixel depending on Ka). Excitation of the dye-doped liquid and imaging of the emitted fluorescence is performed from underneath (unwetted side), in order to avoid spatially and temporally non-uniform beam stirring and lensing in the first instance, and to limit image distortions in the latter. The imaging planes of the two cameras are mapped using a calibration graticule target immersed inside the liquid solution. Optical distortion corrections are performed using the same arrangement, and a pinhole model available in LaVision DaVis.

PLIF/PTV data acquisition and processing
A sample PLIF frame after correction for perspective distortion is depicted in figure 3(a). Near the solid-liquid interface, reflections from the substrate blur the PLIF signal locally, while the fluorescence emitted by the illuminated liquid volume is reflected about the liquid-gas interface and back to the camera. The location of the solid-liquid interface is obtained using an edge-detection algorithm, while for the gas-liquid interface, intercepts between linear fits to maximum signal gradients and reflection intensity profiles are employed as estimates of the liquid boundaries. The processed profiles are then used to mask out the reflection regions and to produce binarised images (see figure 3b), where a signal intensity of unity corresponds to the film domain and a signal intensity of zero is ascribed to the rest of the image. The binarised images are used to mask out reflections from the glass surface underneath the liquid domain and about the gas-liquid interface in the raw particle images (see figure 3c,d). A four-pass cross-correlation algorithm is implemented on the masked-particle-image sequences in order to generate 2D velocity-vector maps (PIV calculation). Following this step, individual particles are tracked (PTV calculation) by employment of the obtained PIV results as reference estimators of the velocity field. There are a number of reasons for employing PTV rather than PIV in this study. To begin with, PTV offers an eight-fold increase in the spatial resolution of the velocity measurement based on the selected velocity-vector calculation settings (which are determined based on the optical set-up, the size of imaged region of the flow and the employed particles). A second reason is that PIV suffers from large bias errors in the presence of strong velocity gradients, as the motion of particle groups is averaged within the PIV interrogation window. Along a smooth film region and near the gas-liquid interface this effect is of little consequence, and PTV and PIV produce similar results. Near the wall, however, as well as between the wave crests and troughs, and along the capillary-wave regions, where the velocity varies strongly over a limited spatial domain, PIV falls short of providing credible quantitative results. It should finally be noted that the presence of particle-image reflections within the masked liquid domain constitutes a commonly encountered error source affecting the PTV velocity data (Kaehler,  Scharnowski & Cierpka 2012). In our experiments, a spatial extent corresponding to only approximately 5 %-7 % of the examined films is strongly affected by such artefacts, and consequently, the impact on bulk velocity measurements is minimal. For a detailed examination of the velocity distributions underneath the waves, the practice of averaging PTV maps corresponding to the same spatial region is adopted. The procedure is initialised by selecting a reference film-thickness profile pertaining to the desired topology, and cross-correlating it with all thickness traces from the same data set. Signal pairs satisfying a maximum displacement condition are then repositioned and averaged (a typical result is shown in figure 4a). The same translation values are used to rearrange and average out the corresponding instantaneous PTV velocity maps, producing averaged velocity fields such as the one displayed in figure 4(b).

Validation
Independent experiments are conducted in order to assess the validity and accuracy of the combined optical methodology. Film thicknesses from flat films (Ka = 14) are compared to predictions from the Nusselt flat-film solution and direct micrometrestage measurements. In both cases the mean relative-deviation amounts to below 1 %. The same uncertainty value is believed to apply to wave-locked average measurements of mean film thickness; however, thinner films ensue in that case, and consequently a worst-case uncertainty of 2 % is a more representative estimate. An estimate of the film-thickness measurement noise is obtained by calculating standard deviations over a wide range of flat-film datasets. This amounts to <30 µm, and therefore, a mean relative error of 2.5 % is quoted for the instantaneous film-thickness measurements.
Relative deviations are also calculated between PTV-derived interfacial and bulk velocities, and corresponding analytical results of the Nusselt flat-film solution, with average relative deviations corresponding to 3.2 % for both velocities. Standard deviations of local velocity in flat films (conservative estimates of PTV measurement noise) are observed to decrease with increasing distance from the wall; at the gas-liquid interface, they typically amount to less than 1 % of the local mean value, while halfway across the film thickness they amount to approximately 5 %-6 % of the mean value.
Solitary waves on falling liquid films in the inertia-dominated regime 501 3.2. DNS methodology The applied DNS methodology is based on a general-purpose finite-volume framework for interfacial flows (Denner & van Wachem 2014a), which resolves the dynamics of the full two-phase system, i.e. the liquid film, the gas-liquid interface and the gas phase. The flow, which is assumed to be incompressible and isothermal, is governed by the momentum equations and the continuity equation where t represents time, u is the velocity, p is the pressure, ρ is the density, µ is the dynamic viscosity, g is the gravitational acceleration and f σ represents the contribution of the surface tension.

Numerical framework
The primitive variables are obtained using a coupled, implicit finite-volume framework with collocated variable arrangement (Denner & van Wachem 2014a). The momentum equations (3.2) are discretised using a second-order backward Euler scheme for the transient term, while the convection, diffusion and pressure terms are discretised using central differencing. The continuity equation (3.3) is discretised using a balanced-force implementation of the momentum-weighted interpolation method, as proposed by Denner & van Wachem (2014a), which couples pressure and velocity.
The volume-of-fluid (VOF) method (Hirt & Nichols 1981) is adopted to capture the gas-liquid interface. In the VOF method the local volume fraction in each cell is represented by the colour function γ , defined as γ = 0 in fluid a and γ = 1 in fluid b. Hence, the interface is located in mesh cells with a colour function value of 0 < γ < 1. The colour function is advected by the linear advection equation which is discretised using a compressive VOF methodology (Denner & van Wachem 2014b). Assuming surface tension is a volume force acting in the interface region, the surface force per unit volume is described by the continuum surface force (CSF) model (Brackbill, Kothe & Zemach 1992)  The properties of the liquid film are directly taken from the experimental measurements (see § 3.4) and the gas phase is assumed to be air with density ρ g = 1.205 kg m −3 and viscosity µ g = 1.82 × 10 −5 Pa s. At the bottom (liquid-side) wall a no-slip condition is enforced and at the top (gasside) boundary a free-slip wall is applied. A monochromatic forcing is imposed at the domain inlet by periodically changing the mass flow at a given frequency f and amplitude A from the mean. A semiparabolic velocity profile is prescribed for the liquid phase at the inlet and a spatially invariant velocity is prescribed at the inlet of the gas phase, given as The co-flow in the gas phase is applied to prevent backflow at the domain outlet, which can significantly affect the computational performance as well as the accuracy of the results, but has a negligible effect on the evolution of the film flow, since the momentum of the gas phase is approximately three orders of magnitude smaller than the momentum of the liquid film. A constant film height h(x = 0) = h N is prescribed at the inlet. An open boundary condition is applied at the domain outlet, as previously used for instance by Nosoko & Miyara (2004) and Demekhin et al. (2007) (Kalliadasis et al. (2012, appendix F3) gives details on how to apply an outlet boundary condition).

LD modelling
We adopt the LD formulation corresponding to the two-field second-order model derived by Ruyer-Quil & Manneville (1998), a set of two coupled nonlinear equations for the evolution of the interface thickness h(x, t) and local flow rate q(x, t): where the Shkadov scaling (Shkadov 1977) has been used, giving rise to the reduced Reynolds number δ defined in (2.7) and the viscous dispersion number η defined in (2.8). The above model is obtained by combining the long-wave expansion up to second order with a weighted residual technique based on a Galerkin projection in which the velocity field is expanded onto a basis with polynomial test functions (Ruyer-Quil & Manneville 1998, 2000 (this formulation was extended to 3D in the study by Scheid et al. (2006)). It is worth remarking that the model contains the second-order viscous dispersion terms (those gathered in the second line of (3.7)), which are controlled by η, often neglected in thin-film studies. These terms affect both the amplitude and frequency of the capillary ripples at the front of solitary waves (Ruyer-Quil & Manneville 2000). The first-order model, obtained by setting η = 0 in the second-order one already corrects the shortcomings of the first-order model obtained by Shkadov (1967), with the principal one being an erroneous prediction of the instability onset, while the second-order model is in good agreement with Orr-Sommerfeld for a sufficiently large region of Re. Also, as was demonstrated by Pradas et al. (2011), who developed a coherentstructure theory for falling films, the second-order terms play an important role in the interactions between solitary pulses precisely because they affect both the amplitude and frequency of the front-running capillary ripples, and hence they are crucial for an accurate description of pulse interaction. So their influence is in fact linear, but they can have profound consequences on the nonlinear dynamics of the film and the wave-selection process in the spatiotemporal evolution. This is because solitary waves interact through their tails which overlap, i.e., the capillary ripples and their amplitude and frequency will affect the separation distance between the pulses: for instance, smaller-amplitude ripples will allow for more overlaps between the tails of neighbouring solitary waves, thus decreasing their separation distance. This in turn will affect the average separation distance between the waves when the system reaches its permanent quasiturbulent regime, and hence the density of the solitary waves.
Equations (  The corresponding Kapitza numbers are Ka = 14 for cases A, Ka = 85 for cases B and Ka = 346 for cases C. In all cases the inclination angle of the substrate is β = 20 • with respect to the horizontal. The film is perturbed with a monochromatic, periodic change of the flow rate, with frequency f and amplitude A, applied at the inlet of the test domain. The applied forcing frequencies are well above the critical frequency, as reported by Argyriadi, Serifi & Bontozoglou (2004), below which parasitic pulses prevent the development of a regular wave train. Based on the classification of Ooshida (1999), all considered cases are situated in the drag-inertia regime, δ 1.
To allow a direct comparison of the experimental and numerical results, all cases are analysed at downstream distance x m = 0.256 m from the inlet of the test domain, as dictated by the employed experimental arrangement. The experiments are conducted first by choosing the amplitude of the monochromatic inlet forcing so that fully developed, quasistationary interface waves are obtained at the measurement position x m = 0.256 m. Subsequently, the numerical simulations are conducted using the measured fluid properties and flow rate as well as the imposed monochromatic forcing frequency and forcing amplitude in the experiments, thus establishing a direct comparison under the same conditions. The presented velocity profiles are taken at the positions illustrated in figure 6, where d is the distance between the crest of the solitary wave and its preceding trough. Figures 7 and 8 show the spatial development of the film height and the streamwise velocity at the interface for cases A2 and C6, respectively, obtained with DNS. The velocity, similar to previously reported comparisons between experiments and DNS by Dietze & Ruyer-Quil (2013), using the open-source software GERRIS, and by Doro & Aidun (2013) and Rohlfs, Pischke & Scheid (2017), using the open-source software OPENFOAM. The LD model is also in very good agreement with both experimental and DNS results provided that δ 12, but fails to capture the amplitude and frequency of the front-running capillary ripples for δ > 12 (which in turn affect  figure 11 for all considered cases with forcing frequency f = 7 s −1 , are generally in satisfactory agreement with each other. In particular, the same trends can be extracted with respect to the change in film height and wave speed as Re increases. An increase in Re causes initially an increase in wave height, irrespective of the Kapitza number Ka, both with respect to the Nusselt flat-film height h N and the film height in the preceding trough h min . Eventually, the normalised maximum film height shown in figure 11(a) reaches a maximum with respect to Re. Increasing Re further ensues a decrease in maximum film height, as already noted in previous studies (Chakraborty et al. 2014;Denner et al. 2016). Interestingly, the observed trend persists for the wide range of Kapitza numbers Ka considered in this study, with Ka = 14-346 in figure 11, both in the experimental measurements as well as the DNS results. The minimum film height h min shows a similar yet inverse behaviour compared to the maximum film height as far as the influence of Re is concerned, similar to the observations reported by Denner et al. (2016) for falling water films. Figure 12 compares the profiles of the dimensionless streamwise velocity u/u N in the liquid film at different positions in the solitary wave, for selected cases. The agreement between experimental measurements and DNS results is in general very good, especially for the small Reynolds numbers, cases A1 (Re = 5.1) and B1 (Re = 7.5). At larger Reynolds numbers the flow field is strongly dependent on the local wave shape and can vary rapidly in space, as shown by the experimental measurements of Charogiannis et al. (2017) for instance. Hence, in particular at the front of the solitary wave (see figure 12c), small differences in local film height and measurement position can lead to noticeable differences in the velocity profile. Note that the velocity profile is parabolic in the tail and underneath the crest of the solitary wave, even for the case with the highest Reynolds number (case C6, Re = 77). This agrees with previously reported numerical results (Malamataris et al. 2002;Gao et al. 2003;Malamataris & Balakotaiah 2008) and experimental measurements (Adomeit & Renz 2000;Moran et al. 2002;Charogiannis et al. 2015). Interestingly, the self-similarity of the velocity profile in the laboratory frame of reference is unaffected by the onset of a recirculation region in the moving frame when the maximum flow velocity exceeds the wave speed (for instance, in case C6). Moving further downstream to the front of the solitary wave, position X3, the velocity profile increasingly departs from the parabolic one that is observed for the tail and the crest of the solitary wave, as clearly seen in figure 12(c). The limitations of the applied PIV/PTV method to measure the velocity fields in our experimental arrangement prohibit a direct comparison of the velocity profiles obtained experimentally with those from DNS in the trough of the solitary wave (position X4). In summary, carefully conducted experiments and DNS using state-of-the-art methods show good agreement when predicting harmonically excited solitary waves Solitary waves on falling liquid films in the inertia-dominated regime 509 on falling liquid films in the inertia-dominated regime and not only for flows with small Reynolds numbers reported in previous studies (for example, Gao et al. (2003), Dietze & Ruyer-Quil (2013)). Note that this comparison of experiments and DNS is more detailed than previously published comparisons (notably Dietze & Ruyer-Quil (2013), Doro & Aidun (2013), Rohlfs et al. (2017)), as it comprises the film height, film velocity at the interface, the velocity profile in the liquid film at different positions, and the wave speed, across a range of Reynolds and Kapitza numbers. This allows one to draw conclusions based collectively on both the experimental measurements and the DNS results. And although the LD model is applied outside what is usually considered to be its range of validity (δ 10), its predictive accuracy with respect to the maximum and minimum film height, and consequently the wave height, is impressive and somewhat unexpected. This might be explained by the predominantly parabolic velocity profile in the tail and underneath the crest of the solitary wave, which accounts for the largest part of the solitary wave, even after the onset of flow recirculation. Re. Thus, the wave height should increase with increasing Reynolds number, instead of stagnating or even reducing, as suggested by the presented numerical results and experimental measurements. Surface tension stabilises the wave, yet the absolute value of the interface curvature of the main wave hump is decreasing for increasing Reynolds number once the wave height does not further increase, since the wavelength of the solitary wave is λ ∼ √

Mechanism for stabilising the wave height
Re sin β (Denner et al. 2016). Hence, the surface tension acting on the solitary wave decreases with increasing Re. Consequently, there must be an alternative mechanism acting to stabilise the height of the solitary wave against the destabilising effect of the increased inertia. Figure 13 shows the maximum film height h max , normalised by the Nusselt height h N , alongside the corresponding maximum flow rate q max , normalised by the Nusselt flow rate q N , for cases with Ka = 346 (cases C) and f = 10 s −1 , in the Reynolds number range Re = 10-85, obtained with DNS. The instantaneous local flow rate is defined as where h(x) and u(x) are the local film height and flow velocity, respectively. Both h max /h N and q max /q N exhibit a clear maximum with respect to the Reynolds number at Re c ≈ 50, and decrease for further increasing Reynolds number Re > Re c , as is evident in figure 13. A similar behaviour is observed comparing h max /h N and q max /q N of a water film falling on a vertical substrate, as previously considered by Denner et al. linear relationship between local flow rate and local film height for a given solitary wave, with a slope that corresponds to the wave speed. Hence, the simultaneous increase and decrease of h max /h N and q max /q N observed in figures 13 and 14 is to be expected. Moreover, the shape of the solitary wave is governed by the balance of inertia, represented by the dynamic pressure p dyn = ρu 2 /2, and surface tension, represented by the capillary pressure p σ = σ /h N (Denner et al. 2016). The relevant characteristic parameter is, thus, the Weber number We = 2p dyn /p σ (see (2.4)), which can be rewritten with respect to position x as based on the local flow rate given in (4.1). Since the density ρ l and surface tension coefficient σ are assumed to be constant, a change in flow rate leads to a change in film height at a given Weber number. Note that the dynamic pressure p dyn represents the (streamwise) translational kinetic energy E u = mu 2 /2 per unit volume of the liquid film, where m is the fluid mass. A plot of the maximum film height against the maximum flow rate, presented in figure 15(a), confirms the correlation between h max /h N and q max /q N established above, but it also shows a reduction of the maximum flow rate corresponding to a given maximum film height for Re > Re c . This difference between q max at a given h max for Re < Re c and the corresponding q max for Re > Re c suggests a reduction of the effective inertia acting on the solitary wave for Re > Re c , with Re c ≈ 50 for the conditions analysed in figure 15(a). Interestingly, this reduction of effective inertia coincides with the development of flow recirculation in the main wave hump with respect to the moving frame of reference when the maximum flow velocity surpasses the wave speed, u max > c, as observed in figure 15(b). The flow recirculation in the moving frame of reference increases the angular velocity Ω and, as a result, increases the angular kinetic energy E Ω = IΩ 2 /2, where I is the moment of inertia. Due to the conservation of energy, the increase of angular kinetic energy E Ω reduces the translational kinetic energy E u . This transfer of energy leads to the observed reduction in effective inertia of the film flow. It is worth noting that the increase in Ω corresponds to an increase in vorticity ω = 2Ω (assuming incompressible flow and neglecting high-order contributions). However, an analysis of the vorticity distribution in the solitary waves of the considered cases has not provided conclusive evidence that the increase in viscous dissipation plays a significant part in the observed reduction of the effective inertia acting on the solitary wave. Note that for the flow conditions studied in figure 15, a theoretical criterion derived by Rohlfs & Scheid (2015) predicts the onset of flow recirculation for Re circ ≈ 20, which is in reasonable agreement with the observed onset of flow recirculation at Re circ ≈ 25. A possible explanation for this difference in Re circ is the reduced bulk velocity compared to the Nusselt solution for increasing wave height reported by Charogiannis et al. (2017), that shifts the critical local flow rate for the onset of flow recirculation to a higher Reynolds number. Aside from the reduction of (streamwise) translational kinetic energy, the onset of flow recirculation in the moving frame of reference also leads to a downward motion (in the wall direction) of the flow in the main wave hump. Before the onset of flow recirculation, the propagation of the solitary wave with c > u max leads to an upward motion (away from the wall) downstream of the wave crest, as observed in figure 16(a). However, after the onset of flow recirculation (c < u max ), the flow moves downwards under the crest and upwards in the tail of the solitary wave, as observed in figure 16(b,c) and as a consequence of the flow recirculation. This change in cross-stream motion contributes to the observed saturation of the film height after the onset of flow recirculation. In general, this effect is similar to the long-wave instability mechanism that generates the solitary waves; at the onset of the instability, the phase difference between vorticity and the instability wave is destabilising if c > u max , yet it is stabilising if c < u max (Kelly et al. 1989). This mechanism suggests that, in general, the fluid motion in the bulk phase has a stabilising effect when c < u max (i.e. flow recirculation in the moving frame of reference), as confirmed by the presented results.
In summary, the presented results show that the flow recirculation in the moving frame of reference is responsible for the stabilisation, or even reduction, of the height Solitary waves on falling liquid films in the inertia-dominated regime 513 of the solitary waves at large Reynolds numbers. The flow rate and the maximum velocity in the film exhibit their highest values at (approximately) the same Reynolds number at which the film height reaches its maximum. With respect to the link between film height and flow rate, this is to be expected based on the recent study of Charogiannis et al. (2017), reporting a linear relationship between film height and flow rate. However, the connection of the film height (and flow rate) with the maximum flow velocity and the recirculation in the moving frame of reference has not been reported before. This finding is of particular interest for the optimisation of the heat and mass transport characteristics of wavy falling liquid films, which has previously been shown to be critically dependent on the characteristics of the solitary waves (Mathie et al. 2013;Markides et al. 2015) and with flow recirculation in the moving frame of reference typically improving the heat and mass transfer within the liquid film (Roberts & Chang 2000;Albert et al. 2014).

Flow reversal
Flow reversal is driven by a strong positive streamwise pressure gradient that results from the rapid change in interface curvature between a solitary wave and its preceding capillary ripple (Dietze et al. 2008(Dietze et al. , 2009 or its preceding solitary wave (Doro & Aidun 2013). Furthermore, previous studies (Dietze et al. 2008(Dietze et al. , 2009Doro & Aidun 2013;Rohlfs & Scheid 2015) found that a larger magnitude of the interface curvature in the trough and at the preceding capillary ripple, as well as a higher inclination angle of the substrate, favour the onset of flow separation and reversal. However, Chakraborty et al. (2014) reported that flow reversal is limited by lower and upper limits with respect to the (reduced) Reynolds number.
The profiles of velocity and streamwise pressure gradient for cases B4 and C2, with a Reynolds number of Re = 26.9 and Re = 28.5, respectively, and a forcing frequency of f = 10 s −1 , are shown in figure 17. The higher inertia relative to surface tension of case B4 (We = 1.98) compared to case C2 (We = 0.53) means the high pressure gradient is concentrated near the interface in case B4, whereas in case C2 the region of high pressure gradient extends to the wall. As a result, flow reversal is observed in case B4 but not in case C2. This happens although the curvature of the trough in cases B4 and C2 is comparable, with κ min = −860 m −1 in case B4 and κ min = −878 m −1 in case C2, and the difference in surface tension between the two cases is merely 2 %. Hence, the interface curvature can be ruled out as the source of the observed differences, indicating that the larger influence of inertia compared to surface tension of case B4 (higher We) is the dominating factor. A similar behaviour can be observed in figure 17 for cases B2 (We = 0.49) and C1 (We = 0.13), which have a Reynolds number of Re = 11.7 and Re = 12.4, respectively, where flow reversal is observed in case C1 but not in case B2.
Interestingly, the magnitude of the flow reversal region as well as the size of the separation eddy are considerably larger in case C2 than in case C1, although the forcing frequency and the Froude number Fr = √ 3Re/ cot β, which is Fr = 3.7 for case C1 and Fr = 5.6 for case C2, mean flow separation should be less likely in case C2, based on findings of previous studies (Chakraborty et al. 2014;Rohlfs & Scheid 2015). The front of the main wave hump in case C2 (d/h N = 4.79) is significantly steeper than in case C1 (d/h N = 7.55), also seen in figure 18, leading to a considerable increase in negative streamwise pressure gradient immediately upstream of the wave trough, which generally counteracts flow separation. However, the shorter wavelength combined with the higher Reynolds number of case C2 (λ/h N = 39.7, Re = 28.5) compared to case C1 (λ/h N = 52.3, Re = 12.4), leads to larger capillary waves for case C2. As a result, the pressure underneath the first capillary ripple is noticeably higher in case C2 than in case C1, as observed in figure 18, which results in a larger streamwise pressure gradient between the trough preceding the main wave hump and the crest of the first capillary ripple, and, consequently, promotes flow separation and reversal. The effect of the pressure changes in the streamwise direction on the velocity underneath the trough can clearly be seen by the strongly curved and clustered isocontours of the streamwise velocity in figure 18.
In the above cases, the low-pressure region underneath the wave trough cannot extend to the low-velocity region near the wall in cases where inertia is high compared to surface tension, limiting the high streamwise pressure gradients to the close vicinity of the interface. Hence, the substantially reduced streamwise pressure gradient near the wall is not able to initiate flow separation, i.e. the direction of the flow remains unchanged. In addition, the front of the solitary wave becomes steeper as the Reynolds number increases (Denner et al. 2016), which increases the curvature of the wave trough and leads to a larger streamwise pressure gradient underneath the trough that promotes flow reversal. Yet, the steeper front of the solitary wave also increases the pressure difference between the main hump and the wave trough, as observed in figure 18. Thus, flow reversal is eventually arrested when the Reynolds number is sufficiently high, because the (negative) pressure gradient upstream of the wave trough dominates the (positive) pressure gradient downstream of the wave trough. As a consequence, the parameter space in which flow reversal occurs is governed by a non-trivial interplay of a variety of parameters, including inclination angle and frequency, as suggested by previous studies (Chakraborty et al. 2014;Rohlfs & Scheid 2015), as well as wave steepness, the size of capillary ripples and the balance of inertia and surface tension.

Conclusions
Falling liquid films exhibit rich spatiotemporal dynamics governed by a complex interplay between inertia, viscosity, surface tension and gravity. Despite several decades of intensive research, this interplay is not fully understood and important aspects of falling film dynamics, especially for inertia-dominated falling films, have not yet been resolved. In this study we undertook a systematic investigation comprising advanced computations and experimentation aimed at addressing open questions associated with the inertia-capillary interactions that govern the shape and propagation of solitary waves on falling liquid films in the inertia-dominated regime.
The overall agreement between experiments and DNS is very good for all quantities of interest, such as film thickness, velocity profile and interface velocity, and within the expectations based on the experimental tolerances and the numerical assumptions. The applied LD model also performs very well and shows excellent agreement with experiments and DNS for reduced Reynolds numbers δ 12. For higher δ the results of the LD model increasingly deviate from the experimental measurements and DNS results, which can be attributed to the boundary-layer approximation on which the LD model is founded; for instance, LD models are not able to predict the pressure fields shown in figure 18 that also include a significant hydrodynamic (due to the complex flow field) contribution. Nevertheless, the LD model is able to predict accurately the maximum and minimum film heights even for inertia-dominated flows with δ 12, a previously unknown, very useful and interesting feature of LD models. A possible explanation for this behaviour of LD models is the parabolic velocity profile in most parts of the film, in particular under the crest of solitary waves (and even in the inertiadominated regime and after the onset of flow recirculation).
Careful interrogation of the wave characteristics shows a stagnation/reduction of the maximum film height at high Reynolds numbers. While in the present study, as well as in the previous study of Denner et al. (2016), this stagnation/reduction of the maximum film height has been observed for trains of travelling solitary waves with constant forcing frequency, a similar behaviour was reported by Chakraborty et al. (2014) for a single solitary wave travelling on an otherwise flat film (corresponding to an infinite wavelength). The presented results reveal that this reduction in maximum film height coincides with a reduction in the maximum flow rate and suggest that a reduction of the effective inertia acting on the solitary wave hump is responsible for the stagnation/reduction of the height of solitary waves at high Reynolds numbers.

F. Denner and others
This reduction of effective inertia is found to be the result of flow recirculation in the main wave hump in the moving frame of reference. The observed maximum in flow rate at a critical Reynolds number, which coincides with the maximum film height, is an important finding as far as exploitation of falling liquid films for heat and mass transfer applications is concerned, since previous studies (Roberts & Chang 2000;Mathie et al. 2013;Albert et al. 2014;Markides et al. 2015) suggested a significant impact of the film height, the flow rate and the recirculation on the heat and mass transfer characteristics of film flows. Regardless of the Reynolds number and regardless of whether or not flow recirculation in the moving frame is present, the velocity profile in the tail of the solitary waves and underneath the wave crest is found in both the experiments and DNS to have a parabolic, self-similar shape.
At the trough preceding the solitary wave, the velocity profile is strongly influenced by the low-pressure region underneath the trough, due to the convex shape of the interface and the resulting positive pressure gradient downstream of the trough. If the conditions are favourable, the flow underneath the trough preceding the solitary wave separates from the wall and flow reversal ensues. However, for sufficiently high inertia, the streamwise pressure gradient from the convex interface shape cannot reach the low-velocity region near the wall; thus, the magnitude of the streamwise pressure gradient is insufficient to initiate flow separation and flow reversal. This agrees with the observation of an upper limit for the flow reversal with respect to the Reynolds number reported, albeit without studying the underlying mechanism in detail, by Chakraborty et al. (2014). In addition to inertia, our results also confirm that the pressure resulting from the curvature of the first capillary ripple is a dominant factor for the flow separation observed in front of the solitary wave, as initially proposed by Dietze et al. (2008). The change in flow direction associated with the flow reversal occurring underneath the wave trough may enhance the heat and mass transport characteristics and, hence, understanding the reasons underpinning the upper (and lower) limit of flow reversal with respect to the Reynolds number is of direct practical relevance.