Georg F Dietze

G Lavalle

C Ruyer-Quil

Falling liquid films in narrow tubes: occlusion scenarios

Keywords:

We study a gravity-driven wavy liquid film falling down the inner surface of a narrow cylindrical tube in the presence of an active core gas flow. We employ the model of Dietze and Ruyer-Quil (

Introduction

We consider the configuration in figure 1, a thin liquid film falling down the inner surface of a vertical cylindrical tube of radius R ⋆ under the action of the gravitational acceleration g. The film of liquid (denoted with the subscript l) is in contact with a laminar flow of gas in the core (subscript g). Both fluids are assumed to be Newtonian with constant densities ρ l and ρ g , dynamic viscosities µ l and µ g , and surface tension σ. Also, we assume the arrangement to be axisymmetric. We denote the film thickness h, the core radius d, its spatial average d, and the flow rates q l and q g . Their dimensional † Email address for correspondence: dietze@fast.u-psud.fr (low-viscosity silicone oil), Re=14.6, Reg=-17.27 (limit point marked by asterisk in panel 13a).

Streamlines in the wall-fixed reference frame within the liquid (blue lines) and gas (red lines).

counterparts will be distinguished by a star, e.g. d ⋆ , and their counterparts in the corresponding flat-film primary flow will be denoted with a subscript zero, e.g. d 0 .

We focus on conditions, where: (i) the effect of gravity, quantified through the Bond number Bo=ρ l g R ⋆2 /σ, is at least comparable to that of capillarity (Bo 1); (ii) the liquid film's inertia, quantified through the Reynolds number Re=q ⋆ l0 ρ l /(πR ⋆ µ l ), is not necessarily negligible (Re 20); and (iii) the effect of the gas flow can be relevant, but its Reynolds number Re g =q ⋆ g0 ρ g /(πR ⋆ µ g ) remains moderate (|Re g |<20).

We investigate this flow with the weighted residual integral boundary layer (WRIBL) model of [START_REF] Dietze | Films in narrow tubes[END_REF], which was previously applied only to gravity-free films. The model accounts for gravity, inertia, full interfacial curvature, axial viscous diffusion, and full inter-phase coupling. Our model computations are confronted with our own direct numerical simulations (DNS) in appendix A and panel 5d.

We focus on the occlusion of the narrow tube due to the interfacial instability of the falling liquid film, as observed experimentally by [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] and [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF]. Occlusion is relevant for chemical engineering applications such as tubular falling film micro-reactors [START_REF] Seebauer | Tube bundle falling film microreactor for performing gas liquid reactions[END_REF]. It also occurs in the human pulmonary airways due to the collapse of the mucus-serous film lining their inner surface [START_REF] Grotberg | Respiratory fluid mechanics[END_REF]). In the lower generations of bronchioles, the effect of gravity and possibly inertia is not necessarily negligible [START_REF] Kamm | Is airway closure caused by a liquid film instability?[END_REF].

The falling liquid film is subject to two instability mechanisms causing interfacial deformations: the Plateau-Rayleigh mechanism [START_REF] Plateau | Recherches expérimentales et théorique sur les figures d'équilibre d'une masse liquide sans pesanteur[END_REF][START_REF] Rayleigh | On the instability of cylindrical fluid surfaces[END_REF][START_REF] Goren | The instability of an annular thread of fluid[END_REF], due to the azimuthal curvature of the film surface; and the Kapitza mechanism [START_REF] Kapitza | Wave flow of thin layer of viscous fluid (in Russian)[END_REF]), due to inertia, which sets in as a result of the gravity-driven mean flow [START_REF] Brooke Benjamin | Wave formation in laminar flow down an inclined plane[END_REF][START_REF] Yih | Stability of liquid flow down an inclined plane[END_REF]. The strength of this mean flow determines whether the absolute nature of the Plateau-Rayleigh instability or the convective nature of the Kapitza instability dominates [START_REF] Duprat | Absolute and convective instabilities of a viscous film flowing down a vertical fiber[END_REF]. For the liquids studied here, we find that occlusion is dictated by absolute instability when Bo ≪ 1, whereas convective instability is responsible for occlusion when Bo 1.

In the case of convective instability, the liquid film can produce travelling surface waves, which remain stationary in their reference frame. Gravity plays an important role in shaping these waves, compressing the leading and elongating the trailing wave front, which results in asymmetric tear-shaped wave humps [START_REF] Dietze | On the Kapitza instability and the generation of capillary waves[END_REF]. The resulting distortion of the film surface favours variations in axial (stabilizing) over azimuthal (destabilizing) surface curvature. In related systems, this mechanism has been found to saturate the Plateau-Rayleigh instability in the weakly-nonlinear regime [START_REF] Frenkel | Annular flows can keep unstable films from breakup: nonlinear saturation of capillary instability[END_REF], to the extent that the film looks stable to the naked eye [START_REF] Quéré | Thin films flowing on vertical fibers[END_REF]). In our current problem, we find that it greatly extends the region of existence of stronglynonlinear travelling waves and considerably delays occlusion. [START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF] applied the integral boundary layer approach of [START_REF] Shkadov | Wave flow regimes of a thin layer of viscous fluid subject to gravity[END_REF] to model an annular liquid film falling in a vertical cylindrical tube. The author found that the Plateau-Rayleigh mechanism systematically increases the amplitude of travelling waves versus a planar falling liquid film, where only the Kapitza mechanism is active.

He also observed catastrophic growth, when the amplitude of travelling-wave solutions diverged as a function of 1/R. In those instances, the effect of gravity was negligible, i.e.

Bo ≪ 1. By contrast, we focus on occlusion in the limit Bo 1. Moreover, the model of [START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF] did not account for axial viscous diffusion, which we find to greatly precipitate occlusion in high-viscosity films. [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] performed falling-film experiments in long narrow tubes and determined the occlusion onset in terms of Re for a wide range of liquids. However, no information was available on the wave dynamics responsible for occlusion. By reproducing two of these experiments (runs 13 and 20 there) with our own spatio-temporal computations, we provide this missing information, leading us to distinguish two wave-induced occlusion scenarios. In scenario I, waves causing occlusion emerge directly from linear wave selection. In scenario II, they result from the secondary instability of a regular train of travelling waves, triggering a cascade of coalescence/absorption events that produce increasingly long, more dangerous, waves. The secondary instability [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF] and the ensuing cascaded coarsening dynamics (Chang et al. 1996b) are well known phenomena in planar falling liquid films. [START_REF] Jensen | Draining collars and lenses in liquid-lined vertical tubes[END_REF] investigated annular collars travelling on a gravity-driven liquid film lining the inner surface of a cylindrical tube. The collars were constructed from unduloids, i.e. symmetric constant-curvature equilibrium shapes [START_REF] Delaunay | Sur la surface de révolution dont la courbure moyenne est constante[END_REF]. The author studied under what conditions collars may grow to occlude the tube. Our study extends this work in that the collar shape is not forced to be symmetric, but follows implicitly from the evolution equations. This allows for the gravity-induced distortion of collars, which we find can greatly delay occlusion. [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] performed falling-film occlusion experiments in a transparent tube, using a high-viscosity silicone oil. By confronting these with linear stability calculations, the authors concluded that occlusion is caused by spatially growing nonlinear waves and not by absolute linear instability (Bo=11.06 in their experiments). Further, with the help of a first-order long-wave model obtained from asymptotic expansion [START_REF] Benney | Long waves on liquid films[END_REF], the authors constructed travelling-wave solutions. By continuing these at fixed wave length in terms of Re, they identified a limit point (LP), and conjectured that this may signal the occlusion onset in a real system. However, the authors cautioned that further work was needed to verify this w.r.t. their experiments. This work was initiated by [START_REF] Camassa | Traveling waves for a model of gravity-driven film flows in cylindrical domains[END_REF], but their model did not account for axial viscous diffusion, which we show to strongly affect travelling wave solutions, and the wavelength assumed in their continuation was too long compared to the experiment.

In our current study, we have checked the conjecture of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] based on travelling-wave and spatio-temporal computations with our model. We have successfully confronted these computations with the experiments of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] and find that accounting for axial viscous diffusion greatly improves agreement. Based on our computations, we find that the surface waves causing occlusion indeed systematically lie beyond the limit of travelling-wave solutions (figure 4), confirming the conjecture of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF]. This limit is highly sensitive to the considered wavelength, long waves being more dangerous than short ones, as shown by [START_REF] Camassa | Traveling waves for a model of gravity-driven film flows in cylindrical domains[END_REF] and [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF]. Thus, any predictive criterion for occlusion based on travelling-wave solutions must account for the type of waves that actually occur in a real system.

For this, we introduce the upper conservative occlusion bound Re max , which corresponds to the limit point of travelling-wave solutions at the spatially most amplified frequency of linear waves f max . For Re>Re max , the most amplified surface waves, which typically emerge in an experiment, do not possess travelling states and occlusion is certain to occur in a naturally evolving film. We find that the occlusion experiments of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] and experimental run 20 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] correspond to this regime. Occlusion in this case occurs through scenario I. To determine Re max , we have computed the most amplified wave frequency f max in our numerical continuations by solving the linear spatial stability problem along with our nonlinear model equations.

This had not been attempted by [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] and [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF].

At the other end, occlusion is delimited by a lower conservative bound Re=Re 0 , below which travelling wave solutions always exist, no matter how great the wavelength, and occlusion is impossible. The existence of such a limit was discovered by [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF].

We find that the conservative occlusion bounds Re 0 and Re max delimit a region of conditional occlusion, Re 0 <Re<Re max , where occlusion is theoretically possible (for long enough waves) but does not necessarily occur in a real system. Whether it does, depends on the specific wave dynamics that unfolds over the spatio-temporal evolution of the film, and whether the tube is long enough to accommodate this. In this regime, we find that occlusion is caused by scenario II and that it applies to experimental run 13 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF]. It also opens the possibility of preventing occlusion through coherent inlet forcing, i.e. by forcing waves of sufficiently high frequency. We validate this idea based on spatio-temporal computations.

Unless Bo is very large, we find that the height of travelling waves at the limit point Re=Re max is far from reaching the tube radius. In contrast to what [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] conjectured, travelling-wave solutions are thus lost abruptly and not due to the wave height reaching the tube radius continuously. [START_REF] Zhou | Viscoelastic liquid film flowing down a flexible tube[END_REF] investigated viscoelastic liquid films falling in a deformable narrow tube. The authors performed calculations with a model based on the same approach as [START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF], which does not account for axial viscous diffusion. Interestingly, it was shown that occlusion can cause a contraction of the flexible tube. Such events are known to occur in the pulmonary airways [START_REF] Grotberg | Respiratory fluid mechanics[END_REF]) and can damage cells within the capillaries [START_REF] Bian | Experimental study of flow fields in an airway closure model[END_REF]. The authors went on to determine the occlusion onset by varying the liquid flow rate in spatio-temporal computations. This onset was defined as the point at which the computation breaks down due to an occlusion event.

However, such computations cannot attain a representative developed state. For example, when starting from a flat-film initial condition, a single unrealistically large tsunami-like wave usually develops in the early stages. This wave is bound to cause occlusion, but it is not representative of a real system. By contrast, we characterize occlusion based on the bounds Re 0 and Re max of travelling-wave solutions, and our spatio-temporal computations have all been continued until reaching a statistically-developed state. To achieve this, we have allowed our computations to continue past occlusion events by numerically limiting the core radius of an occluded region to a small but finite value.

Liu & Ding (2017) studied a vertically falling glycerol film flowing down a porous cylindrical surface using a lubrication equation. The authors found that porosity precipitates both the absolute instability threshold and the occlusion limit. Using the same approach, [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF] studied the effect of the Marangoni instability on a radiallyheated liquid film falling in a vertical cylindrical tube. The authors found that Marangoni stresses promote/delay occlusion when the film is heated/cooled from the wall. However, their modelling approach does not account for inertia nor axial viscous diffusion. We find that this does not allow to accurately predict occlusion in high-viscosity (because of axial viscous diffusion) or low-viscosity (because of inertia) liquids.

The model employed in our current study extends upon the earlier works of [START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF], [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF][START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 2. air-driven flow[END_REF] Re max (by 74% for the high-viscosity silicone oil). When the viscosity is low (low-viscosity silicone oil), inertia cannot be neglected and is found to precipitate Re max by 20%. We also find that these two effects determine whether the occlusion mechanism is dictated by absolute linear instability or the loss of travelling-wave solutions.

Our model also accounts for full inter-phase coupling between the liquid and gas, and this has allowed us to study the effect of a laminar counter-current gas flow on the occlusion bound Re max . For low-viscosity liquids, where the wave-induced pressure variation in the gas becomes relevant in the liquid force balance, we find that Re max is significantly reduced with increasing |Re g | (by up to 25% in our computations). Thereby, the conditions we have studied differ from previous works. [START_REF] Alekseenko | Primary instabilities of liquid film flow sheared by turbulent gas stream[END_REF] experimentally studied the linear stability of a falling liquid film subject to a co-current turbulent gas flow. Their tube radius was large and thus occlusion could not occur. [START_REF] Camassa | Ring waves as a mass transport mechanism in air-driven core-annular flows[END_REF][START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 2. air-driven flow[END_REF] studied co-current upward liquid/gas flows based on a longwave model, where the effect of the (turbulent) gas flow was accounted for by relaxing the inter-phase coupling conditions in the limit of high gas velocities. By contrast, we focus on laminar (counter-current) gas flows, where our model, which relies on the unrelaxed coupling conditions, behaves well.

For the case of a planar falling liquid film sheared by a confined gas flow, many works have investigated the flooding phenomenon, where the liquid film either occludes the channel [START_REF] Vlachos | Visual observations of flooding in narrow rectangular channels[END_REF], surface waves [START_REF] Tseluiko | Nonlinear waves in counter-current gas-liquid film flow[END_REF] or the liquid flow [START_REF] Trifonov | Counter-current gas-liquid wavy film flow between the vertical plates analyzed using the Navier-Stokes equations[END_REF]) reverse direction, or the liquid film disintegrates into drops [START_REF] Zapke | Countercurrent gas-liquid flow in inclined and vertical ducts -I: Flow patterns, pressure drop characteristics and flooding[END_REF].

Finally, several works are related less directly to the studied problem. [START_REF] Kouris | Dynamics of axisymmetric core-annular flow in a straight tube. I. the more viscous fluid in the core, bamboo waves[END_REF] studied liquid/liquid gravity-driven flows through a vertical pipe. [START_REF] Beltrame | Partial and complete wetting in a micro-tube[END_REF] studied the transition between partial and complete wetting in a microtube. [START_REF] Xu | Trapping and displacement of liquid collars and plugs in rough-walled tubes[END_REF] studied migration of liquid films in tubes with rough walls.

Moreover, many works have dealt with pressure-driven core-annular flows [START_REF] Aul | Stability of a thin annular film in pressure-driven, low-reynolds-number flow through a capillary[END_REF][START_REF] Joseph | Core-annular flows[END_REF] or liquid films falling down fibres [START_REF] Kalliadasis | Drop formation during coating of vertical fibres[END_REF][START_REF] Quéré | Fluid coating on a fibre[END_REF]). In the latter case, occlusion cannot occur, unless the gas phase is confined by an additional concentric hollow cylinder (Wray 2013).

Our manuscript is structured as follows. Section 2 introduces the governing equations, and details of our model ( §2.1), direct numerical simulations ( §2.2), and linear stability calculations ( §2.3). Results are presented in §3, where we start by proving that occlusion in real systems results from surface waves that lie beyond the limit of travelling-wave solutions ( §3.1). In §3.2, we compute the upper and lower occlusion bounds Re max and Re 0 , delimiting the regimes of certain, conditional, and impossible occlusion. For one of the studied cases, we validate these results based on our own DNS of travelling-wave solutions. In sections 3.3 to 3.6, we show how Re max is affected by gravity, axial viscous diffusion, inertia, and the core gas flow. In §3.7, we demonstrate occlusion scenarios I (certain occlusion regime) and II (conditional occlusion regime) by reproducing two of the experimental runs from [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] with spatio-temporal computations.

In §3.8, we show that occlusion in the conditional regime can be prevented through

Mathematical description

We consider the flow in panel 1a, which is rotationally symmetric w.r.t. the tube axis, and is governed by the Navier-Stokes and continuity equations for the liquid (subscript k=l, χ l =1) and gas (subscript k=g, χ g =Π µ /Π ρ ):

d t u k = -∂ x p k + Fr -2 + χ k Re -1 r -1 ∂ r (r ∂ r u) + ∂ xx u , d t v k = -∂ r p k + χ k Re -1 ∂ r r -1 ∂ r (r v) + ∂ xx v , (2.1a) ∂ x u k + r -1 ∂ r (v k r) = 0, (2.1b)
where Re=ρ l UL/µ l =q ⋆ l0 /(πR ⋆ )/(µ l /ρ l ) is the Reynolds number and Fr=U/ √ gL is the

Froude number, all variables having been rendered dimensionless with the following reference scales: the length scale L=R ⋆ , which corresponds to the tube radius, the velocity scale U=q ⋆ l0 /π/R ⋆2 , which corresponds to the surface velocity of the liquid, the time scale T =L/U, and the pressure scale P k =ρ k U 2 , which corresponds to the phase-specific dynamical pressure. The star symbol denotes dimensional quantities, and the subscript zero the flat-film primary flow. Thus, q l0 designates the nominal liquid flow rate.

The set of governing equations is completed by the boundary conditions:

u l | r=R = v l | r=R = 0, ∂ r u g | r=0 = v g | r=0 = 0, (2.2)
and the kinematic/dynamic inter-phase coupling conditions at the film surface r=d:

u l = u g , v l = v g = ∂ t d + u k ∂ x d, (2.3a) τ t l = Π µ τ t g , τ t k = ∂ r u k + ∂ x v k -2∂ x d 1 + ∂ x d 2 -1 [∂ x u k -∂ r v k ] ,
(2.3b)

τ n l + Re p l = Π µ τ n g + Π ρ Re p g + ReWe -1 κ, τ n k = 2 1 + ∂ x d 2 -1 ∂ x d ∂ r u k -∂ r v k -∂ x d 2 ∂ x u k + ∂ x d ∂ x v k , (2.3c)
where τ t and τ n designate the tangential and normal interfacial viscous stresses, Π ρ =ρ g /ρ l and Π µ =µ g /µ l are the density and dynamic viscosity ratios, We=ρ l U 2 L/σ is the Weber number, and κ denotes the total curvature of the film surface:

κ = ∂ xx d - 1 d 1 - 1 2 (∂ x d) 2 .
(2.4)

We will refer to computations based on the full equations (2.1) to (2.3) as direct numerical simulations (DNS). These have been performed for validation purposes and results are reported in appendix A and panel 5d. However, most of our computations were carried out with the weighted residual integral boundary layer (WRIBL) model of [START_REF] Dietze | Films in narrow tubes[END_REF], which we introduce next.

Weighted residual integral boundary layer (WRIBL) model

The WRIBL model consists of three coupled partial differential equations for the core radius d, and the flow rates q l and q g (Dietze & Ruyer-Quil 2015):

∂ x q j -ε j 2π d ∂ t d = 0, (2.5a) {S i ∂ t q i + F ij q i ∂ x q j + G ij q i q j ∂ x d} = Fr -2 (1 -Π ρ ) -We -1 ∂ x [κ] + Re -1 (C jl -Π µ C jg )q j +Re -1 J j q j (∂ x d) 2 + K j ∂ x q j ∂ x d + L j q j ∂ xx d + M j ∂ xx q j , (2.5b)
where the subscripts i and j are to be permuted over the phase indicators for the liquid (i,j=l) and gas (i,j=g). This yields two equations from (2.5a), ensuring integral mass conservation, and one equation from (2.5b), which ensures conservation of momentum.

The RHS terms in (2.5b) account for different driving forces, i.e. gravity, capillarity, and viscous drag. Relating these terms directly to one another yields two alternative dimensionless groups that are also useful for characterizing the studied flow regimes: the Bond number Bo=We Fr -2 =ρ l gR ⋆2 /σ, which relates gravity to capillarity, and the capillary number Ca=We/Re=µ l U/σ, which relates viscous drag to capillarity. Following [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], we will also make use of the Kapitza number Ka=σ ρ -1 l g -1/3 µ -4/3 l , and the Laplace number La=σ ρ l R ⋆ /µ 2 l . We point out however that the list Re, Fr, We, Bo, Ca, Ka, and La, contains only three independent dimensionless groups, e.g. Re, Ka, and Bo, which (along with Π µ and Π ρ ) suffice to fully describe the flow. All other groups can be expressed in terms of these three, e.g. Ca=Re(BoKa 3 ) -1/2 .

The model is completed by an equation for the gas pressure gradient

∂ x p g | r=d : 2 Π ρ Re ∂ x p g | r=d = -Re S i ∂ t q i + F ij q i ∂ x q j + G ij q i q j ∂ x d +Re Fr -2 (1 + Π ρ ) -We -1 Re ∂ x [κ] + (C jl + Π µ C jg )q j +J j q j (∂ x d) 2 + K j ∂ x q j ∂ x d + L j q j ∂ xx d + M j ∂ xx q j , (2.5c)
which is used either to evaluate the pressure distribution a posteriori or to impose a gas pressure drop through an integral condition. The coefficients F ij , G ij , C ij , S j , J j , K j , and M j in (2.5b) and their counterparts marked by an underscore in (2.5c) are known functions of d. They can be found in appendix A of [START_REF] Dietze | Films in narrow tubes[END_REF].

We solve the equation system (2.5) numerically using three approaches: (i) transient periodic computations, where the domain length L corresponds to the wavelength Λ and where conditions. At the liquid inlet, x=0, we fix the core radius d=d 0 of the flat-film primary flow. For the liquid flow rate q l , we superimpose a zero-mean temporal perturbation on the nominal value q l0 :

∂ i x q k x=0 = ∂ i x q k x=Λ and ∂ i x d x=0 = ∂ i x d x=Λ for i=0,
q l (x = 0, t) = q l0 [1 + F (t)] ,
(2.6a) where the function F (t) defines the type of inlet perturbation:

F (t) = ǫ 1 sin(2π f t) + ǫ 2 N k=1 sin(2π k ∆f t + ϕ rand ), ∆f = 2 f c /N.
(2.6b)

The first term in (2.6b) constitutes a harmonic perturbation of frequency f and the second one mimics white noise through a series of N =1000 Fourier modes that are shifted by a random phase shift ϕ rand =ϕ rand (k) ∈ [0, 2π] and that span a frequency range of twice the linear cut-off frequency f c (Chang et al. 1996a). All our computations were run with the same ϕ rand (k) number series, which was generated once and for all with the pseudo random number generator RandomReal in Mathematica (2014). The strength of the two terms in (2.6b) is determined through their amplitudes ǫ 1 and ǫ 2 . When ǫ 1 =0, the inlet perturbation consists only of white noise. This setting will be used to simulate the natural, noise-driven, evolution of a wavy film as it would occur in an experiment (sections 3.1 and 3.7). In §3.8, we will use coherent inlet forcing to prevent occlusion. In those computations, we will set ǫ 1 >0, thus adding a monochromatic harmonic perturbation to the inlet noise.

At the outlet, x=L, we have implemented the boundary conditions of Richard et al.

(2016), which ensure that liquid is always sufficiently drained from the domain. We set

d| N +2 = d| N +1 = d| N , and q l | N +1 =q l0 (R-d| N ) 3/2 /(R-d 0 ) 3/2
, where N corresponds to the last grid point within the domain and N + 1, N + 2 to the two downstream ghost points.

Computations were started from the initial condition d(x, t = 0)=d 0 , q l (x, t = 0)=q l0 .

In our codes, the gas flow rate q g is expressed in terms of the total flow rate q tot : q tot (t) = q l (x, t) + q g (x, t), (2.7)

which is spatially invariant. In our computations, we either fix q tot explicitly, which fixes the nominal gas flow rate q g0 , given a nominal liquid flow rate q l0 . Or, q tot is dynamically adjusted to produce a fixed gas pressure drop ∆p g , imposed through an integral condition on (2.5c). We quantify this pressure drop through the normalized pressure gradient:

Ψ = ∆p g /LFr 2 , (2.8)
where L is the domain length. When Ψ=1, the pressure drop exactly balances the weight of the gas column, and we refer to this as the aerostatic situation. This is most realistic for reproducing experiments where the core gas flow is quiescent.

In our transient computations, the wavy liquid film can tend toward occluding the tube.

In an experiment, such an event would form a liquid plug separated by two gas bubbles.

The topological change occurring in this situation cannot be mathematically represented by our model, which implies continuous fluid layers. Thus, we apply a numerical procedure to allow our computations to continue beyond such events. For this, we level the core radius d to an arbitrary value d crit =0.1 whenever/wherever d d crit . As a result, occluded zones are represented by a thin gas filament connecting two bubbles either side of a liquid plug (figure 3c). Although our approach allows to recover some basic features of these zones, it is a crude and invasive approximation of the actual physics. For example, every time d is reset to d crit , a bit of liquid mass is lost and this can cause the unphysical reopening of liquid plugs (supplementary movies M2 and M3). We stress however, that we are not interested in studying the actual behaviour of occluded zones here. Our current paper focusses on the events leading up to occlusion, and thus our numerical limitation of the core radius d is simply a means to continue our computations until a developed state has been reached. Nonetheless, we propose in appendix E an alternative way of representing liquid plugs in the framework of our WRIBL model, based on an additional repulsive force in the momentum equation (2.5b) that stabilizes the film surface at very small d. This approach is inspired by the representation of contact line problems with film models [START_REF] Thiele | Sliding drops in the diffuse interface model coupled to hydrodynamics[END_REF], and confronting it with full-fledged plug models [START_REF] Ubal | Stability of the steadystate displacement of a liquid plug driven by a constant pressure difference along a prewetted capillary tube[END_REF][START_REF] Suresh | The effect of gravity on liquid plug propagation in a two-dimensional channel[END_REF]) is an enticing prospect for future work.

Direct numerical simulation (DNS)

We have validated our model computations with selected DNS based on the full governing equations (2.1)-(2.4). These were performed with two different codes.

Transient simulations were performed with the finite volume solver Gerris [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF], which uses the volume of fluid [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and continuum surface force [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]) methods to represent the two-phase nature of the flow. We have successfully applied this code to falling liquid films in previous works [START_REF] Dietze | Wavy liquid films in interaction with a confined laminar gas flow[END_REF][START_REF] Dietze | Films in narrow tubes[END_REF][START_REF] Dietze | Effect of wall corrugations on scalar transfer to a wavy falling liquid film[END_REF], where further details are given. In the current work, we have used it to perform axisymmetric transient DNS on a domain spanning the wavelength Λ in streamwise direction and the tube radius R in radial direction, applying periodic streamwise boundary conditions. To allow imposing the gas pressure drop ∆p g

(2.8) in this periodic setting, the streamwise momentum equations (2.1) for the liquid (k=l) and gas (k=g) were rewritten in terms of the modified pressure pk =p k +x ∆p k /Λ, yielding the additional source term Γ k :

Γ k = ∆p k Λ , ∆p l = Π ρ ∆p g .
(2.9)

Thanks to this source term, the governing equations can be solved with a periodicity condition on pk , while the actual pressure p k is subject to a pressure drop, allowing to control the gas flow rate. Our transient DNS were started from a flat film initial condition corresponding to the primary flow. All transient DNS were performed for case 4 in table 1, on a grid of square cells with increments ∆ x =∆ r =2 -7 in the bulk of the fluid phases and ∆ x =∆ r =2 -8 around the liquid-gas interface. The time step ∆ t was dynamically adapted so that w ⋆ ∆ t /∆ x 0.8 in each cell, whereby w ⋆ denotes the magnitude of the local dimensional velocity. We have verified that our simulations are grid independent.

Results of our transient DNS are presented in appendix A, where they are used to validate our model computations.

The second DNS code allows to construct travelling-wave solutions of (2.1)-(2.4) based on a pseudo-spectral approach, which we have integrated into the continuation software AUTO07P [START_REF] Doedel | Continuation and bifurcation software for ordinary differential equations[END_REF]. This has allowed us to verify our numerical continuation results

obtained with the WRIBL model (see panel 5d). In these travelling-wave DNS, the gas phase is not accounted for, i.e. Π µ =Π ρ =0, and thus we have only applied them to case 3, where the effect of the gas is negligible. The underlying numerical procedure is detailed in appendix B.

Linear stability analysis

Linearizing the WRIBL model (2.5) around the uniform base state {d, q k }={d 0 , q k0 } and assuming a small-amplitude perturbation {d ′ , q ′ k }: yields the complex dispersion relation: (2.11) where j and k are to be permuted over the phase indicators l and g, ω designates the (complex) wave frequency, α the (complex) wave number, ε ll =ε gg =-1, and ε lg =ε gl =1.

d(x, t) = d 0 + d ′ (x, t) = d 0 + d exp [i (α x -ω t)] q k (x, t) = q k0 + q ′ k (x, t) = q k0 + qk exp [i (α x -ω t)] , (
α i (d) 1E-6 1E-5 1E-4 1E-3 1E-2 1E-1 1E0 1E1 1E-10 1E-8 1E-6 1E-4 1E-2 1E0 1E2 1E4 AI AI ω max Re
DR(ω, α) = 2π d 0 ω {ε jl C jl + ε jg Π µ C jg } + α {∂ d C jl q j0 -Π µ ∂ d C jg q j0 } -2π d 0 α 2 M j ε jl -α 3 L j q j0 -i ReWe α 2 /d 2 0 -α 4 + i Re 2π d 0 ω 2 (S g -Sl) + ω α (F jl q j0 -F jg q j0 ) -α 2 G jk ,
A temporal stability analysis assumes that waves grow self-similarly in their reference frame. This description implies α∈R and ω=ω r +i ω i , where the temporal growth rate ω i is obtained by solving DR(ω, α)=0. In panel 2a, we have plotted the thus obtained dispersion curve ω i (α) of the temporal growth rate for cases 1, 3 and 4 from table 1.

Lines correspond to the model dispersion relation (2.11), and symbols correspond to our own numerical solutions of the Orr-Sommerfeld equations, following the work of [START_REF] Hickox | Instability due to viscosity and density stratification in axisymmetric pipe flow[END_REF]. Agreement between the two data sets is good. Thus, our model is able to accurately represent the linear wave selection.

In a spatially-evolving film, it is more representative to consider the spatial growth rate α i . We thus solve DR(ω, α)=0 for α i (ω), by assuming ω ∈R and α=α r +i α i . In panel 2b, we have plotted the thus obtained growth rate dispersion curves α i (ω) for the three cases from panel 2a. Crosses mark the point of maximal growth rate and we denote ω max the associated most-amplified angular velocity.

In section 3.2, we will characterize occlusion by determining the limit of travellingwave solutions through numerical continuation. This limit is very sensitive to the wave frequency f . Thus we focus on the spatially most-amplified waves, which are most likely to emerge from linear wave selection in an experiment. For this, we impose the linearly most-amplified frequency f max =2π/ω max in our continuation runs.

To determine f max , which changes with the control parameters varied in these runs, our continuation code was augmented to additionally solve for ω max subject to:

DR(ω max , α) = 0, ∂ ω α i | ω=ωmax = 0,
(2.12)

where ∂ ω α i is obtained from ∂ ω DR=0. In a given continuation run, the linear stability problem (2.12) is solved for the primary flow {q l0 , q g0 }={q l , qg }, which is fixed by the mean flow rates ql and qg of the current travelling wave solution.

As shown by [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], the annular falling liquid film can be subject to absolute instability. Thus, we have checked for the absolute instability threshold (AI) in all our continuation runs. In the context of our spatial stability analysis, where ω ∈R, the AI threshold bounds the solutions of (2.12) in terms of the varied control parameter, e.g.

Re, because the maximum growth rate α i (ω max ) diverges there. For example, panel 2c represents dispersion curves α i (ω) for case 3 at different values of Re, increased toward the AI threshold. The latter is marked by an asterisk on the blue curve, which traces the locus of the spatial growth rate maximum α i (ω max ). We point out that ω has been rescaled in panel 2c.

Panel 2d represents how ω max changes with Re for the three cases from panel 2b. These curves allow to discern the AI threshold (marked by asterisks) in terms of Re. For case 4

(dashed curve), we have truncated our solution at Re=100, beyond which our nonlinear model predictions are no longer reliable. The AI threshold lies beyond this point but, as we will see, the occlusion bound Re max lies well within the represented range (panel 6a).

In our current study, we focus on regimes where Bo 1. In this case, the AI threshold is always situated beyond the limit of nonlinear travelling-wave solutions and thus it does not affect the occlusion bound. However, we will demonstrate that neglecting certain physical effects in the mathematical description, such as axial viscous diffusion (panel 9b) or inertia (panel 12a), can change this.

Results and discussion

In subsections 3.1 and 3.2, we reproduce numerically the visualization experiments of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] and demonstrate that the occlusion of a narrow tube by a wavy falling liquid film is caused by surface waves that do not possess a finite-amplitude travelling state. Depending on how these waves emerge in a spatially evolving film, either directly from linear instability (scenario I), or subsequent nonlinear dynamics (scenario II), different occlusion scenarios are possible. In subsection 3.2, we determine the bounds of these two scenarios based on the numerical continuation of travellingwave solutions with our model (2.5). We obtain an upper conservative bound Re max , above which occlusion in a naturally evolving film is certain, and a lower conservative bound Re 0 , below which occlusion is impossible. These bounds delimit a regime of conditional occlusion, Re 0 <Re<Re max , where occlusion is theoretically possible but does not necessarily occur in a real system.

Case

ρ l (kg/m 3 ) µ l (Pas) σ (mN/m) R ⋆ (mm) Ka (-) Π ρ (-) Π µ (-)
1 970 12.9 21.5 5.0 3.3×10 -3 1.2×10 1. Representative parameters for our computations. Liquid properties correspond to a high-viscosity silicone oil (case 1), glycerol-water mixtures concentrated at 99 % and 89 % by mass (cases 2 and 3), and a low-viscosity silicone oil (case 4), while the core fluid is air with µg=1.8•10 -5 Pas and ρg=1.2 kg/m 3 . Case 1 corresponds to figure 3a in [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], and cases 2 and 3 correspond to experimental runs 20 and 13 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF].

The dimensionless groups are Bo=ρ

l g L 2 /σ, Ca=µ l U/σ, Re=ρ l UL/µ l , Ka=σ ρ -1 l g -1/3 µ -4/3 l , La=σ ρ l L/µ 2 l , We=ρ l U 2 L/σ
, and Fr=U/ √ gL, where L=R ⋆ and U=q ⋆ l0 /π/L 2 , q ⋆ l0 designating the dimensional nominal liquid flow rate. Remax is the upper conservative occlusion bound (see e.g. panels 5a and 5b).

In subsection 3.7, we demonstrate occlusion scenarios I and II, which are associated with the regimes of certain (Re>Re max ) and conditional occlusion (Re 0 <Re<Re max ), by reproducing two of the experiments in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] with our own spatiotemporal computations. Before that, in subsections 3.3 to 3.6, we establish how different physical effects, i.e. gravity, axial viscous diffusion, inertia, and the core gas flow, affect the upper conservative bound Re max , which delimits these two regimes. Except for section 3.6, where we will vary the gas Reynolds number Re g , our computations have been performed for an aerostatic pressure drop, i.e. Ψ=1. This is most representative of experiments in a quiescent gas.

Parameters for our computations are varied around the four cases listed in table 1, which represent real liquids and plausible values of the tube radius. Case 1 corresponds to figure 3a in [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], where the working liquid was a high-viscosity silicone oil. Cases 2 and 3 correspond to experimental runs 20 and 13 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], where glycerol-water mixtures concentrated at 99 % and 89 % by mass were used. Finally, case 4 concerns a liquid film of low-viscosity silicone oil. The four cases cover a considerable range of the Laplace number La=σ ρ l R ⋆ /µ 2 l , which relates capillarity to viscous drag, and the Reynolds number, which we have quantified with Re max . We see that inertia, which drives the Kapitza instability, is significant for case 4. For this case, which is the most challenging, we have validated our model computations with transient DNS in appendix A (comparisons with travelling-wave DNS, for case 3, are reported in panel 5d).

3.1. Nature of the surface waves that cause occlusion [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] performed occlusion experiments in an R ⋆ =0.5 mm and L ⋆ =400 mm transparent tube using high-viscosity silicone oil (case 1 in table 1). In their figure 3, the authors produce photographs of the wavy falling liquid film at different Re beyond the

(a) Re=2.3×10 -4 , q ⋆ l0 =4.8×10 -2 cm 3 /s, ǫ2=0.0002; (b) Re=3.2×10 -4 , q ⋆ l0 =6.7×10 -2 cm 3 /s, ǫ2=0.00005; (c) Re=4.5×10 -4 , q ⋆ l0 =9.4×10 -2 cm 3 /s, ǫ2=0.0001.
experimental occlusion threshold, allowing to discern the nature of the surface waves that cause occlusion. We have reproduced three of these experiments with spatio-temporal computations based on our model (2.5) on a long domain using inlet/outlet conditions (see section 2.1). Experimental noise was mimicked through the noisy inlet perturbation

(2.6), where we have set ǫ 1 =0, and tuned ǫ 2 to match the spatial evolution of surface waves in the experiments (when decreasing/increasing ǫ 2 , the region of linear growth of surface waves is stretched/compressed in streamwise direction and this shifts the emergence of nonlinear effects, such as occlusion, downstream/upstream). Figure 3 shows snapshots of these computations, which are to be directly compared to panels 3a, 3b, and 3c in [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF].

Agreement with the experimental snapshots is good. In particular, the wavelength of the four distinct surface waves at the tube outlet in panel 3b is predicted accurately by our computations. Also, the number of gas bubbles (eleven) in panel 3c, which result from successive occlusion events, agrees with the experiment. Panel 3c and the supplementary movie M1 allow to discern the actual mechanics of occlusion. Surface waves in the upstream portion of the tube grow spatially until their crests reach the tube axis, whereupon individual bubbles of the gas core are pinched off. Our model accurately predicts the length scale of the pinch-off process in comparison with the experiment. We stress that panel 3c is a snapshot and that the occlusion point shifts up and down the tube intermittently as a result of the noisy wave spectrum imposed through the inlet condition.

These results show that our model captures accurately the wave dynamics leading up to occlusion, and this is the focus of our current study. However, the physics of the liquid plugs resulting from occlusion is captured only very crudely by the numerical core radius limitation introduced in section 2.1. For example, the shape of the gas bubbles in panel 3c is quite different from that in the experiment and their length evolves as they travel downstream (see also supplementary movie M1). First results reported in appendix E suggest that these shortcomings can be alleviated with our improved WRIBL model. But, the exact shape of liquid plugs can never be fully captured in the context of lubrication theory, as the film surface slope diverges at their edges.

We focus now on panel 3a, where Re is closest to the experimental occlusion threshold. We have also reproduced numerically the experiment in the L ⋆ =1 m tube and a snapshot of our computation is represented in panel 4a. It confirms that occlusion does indeed occur, in the region x 130. Occlusion results from unbounded spatial growth of the prevailing surface waves, and panels 4b to 4d show how the latter emerge from the noisy inlet condition. In these panels, we have represented the single-sided amplitude spectrum P (f j ) of the discrete Fourier transform F (f j ) applied to core radius time series d(t k ) recorded at different streamwise positions x (marked by dashed lines in panel 4a), leading up to the occlusion region:

P (f j ) = (2/N ) F (f j ) ∀ j N/2 + 1, F (f j ) = N k=0 d(t k ) exp(i 2πj/N k), (3.1)
where N =1.8 • 10 6 is the number of (equidistant) sampling points t k with which we have recorded the time series d(t k ) over a duration T =9.1 • 10 4 , and f j =(j/N )/T are the frequencies of the discrete spectrum. At the most downstream location (panel 4d), the spectrum is organized around the spatially most amplified frequency of linear waves f max in the form of a quite narrow peak. This peak lies beyond the limit f /f max =1.12 for the existence of travelling-wave solutions at the considered Reynolds number Re=2.3 × 10 -4 , which is highlighted by vertical blue lines in panels 4b, 4c, and 4d. This limit has been determined by numerical continuation of travelling-wave solutions, which we will discuss in the next subsection (it is marked by an asterisk in panel 5b). We will see that higher frequency (shorter) waves are safe from and lower frequency (longer) waves are prone to occlusion.

By way of numerical continuation, [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] were the first to show that the existence of travelling-wave solutions is bounded. By increasing Re at fixed wavelength, the authors found a limit point beyond which travelling-wave solutions cease to exist.

They conjectured that the loss of such solutions signals occlusion in an actual experiment.

However, the authors cautioned that further work was needed to verify this conjecture.

This work was initiated by [START_REF] Camassa | Traveling waves for a model of gravity-driven film flows in cylindrical domains[END_REF] Figure 4 proves that the occlusion of a narrow tube by a wavy falling liquid film is caused by surface waves that lie beyond the limit of travelling-wave solutions. Thus, occlusion can be characterized based on the existence span of such solutions and we proceed to this in the following subsection.

Travelling waves: regimes of certain, conditional, and impossible occlusion

The frequency of travelling waves determines their amplitude, which is well known for planar falling liquid films [START_REF] Nosoko | Characteristics of twodimensional waves on a falling liquid film[END_REF] and also holds in our cylindrical configuration. At a given value of Re, the smaller the wave frequency, the greater the wave Using our model (2.5), we quantify it for the real liquids studied here (cases 1-4 in table 1), which require a more complicated modelling approach.

Panel 5a represents curves of travelling-wave solutions for case 1 in terms of the minimal core radius d min , as obtained by varying the liquid Reynolds number Re. On each curve, we have fixed the wave frequency f to a certain fraction f /f max of the spatially most amplified frequency f max , which we determine through the approach outlined in section 2.3. All curves display a limit point in terms of Re, beyond which there are no travellingwave solutions. The existence of such a limit point was discovered by Camassa et al.

(2014), who imposed a constant wavelength in their continuation. This wavelength was too long compared to their experiments, whereas we see here that accounting for the correct wavelength/frequency is quite important. Indeed, the limit point of the curves in panel 5a varies greatly with f /f max . The lower the wave frequency, the smaller the limit value of Re, the more danger of occlusion.

We now introduce the upper conservative occlusion bound Re max , which corresponds to the limit point (LP) of the solid line representing the spatially most amplified travelling waves, i.e. f /f max =1. For Re>Re max , the linearly most amplified surface waves, which are most likely to emerge in an experiment, do not possess travelling states and occlusion is certain to occur in a naturally evolving film. The experiment of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], which we have reproduced numerically in figure 4, corresponds to this regime. Occlusion in this case occurs through scenario I, which we will discuss further in section 3.7 based on spatio-temporal computations. bound Re=Re max =1.9×10 -4 and its limit point thus coincides with the maximum of the growth rate dispersion curve α i (f ). We stress that determining Re max requires simultaneously computing the most amplified wave frequency f max in our numerical continuations, as described in section 2.3. This had not been attempted in the previous studies of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF] and [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF].

Reducing Re below Re max moves the limit point of travelling-wave solutions to lower values of f . Thus, the occlusion range is increasingly confined to (very long) waves of very low frequency, until it vanishes completely at Re=Re 0 (red dot-dashed curve in panel 5b). This limit, the existence of which was discovered by [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF], constitutes a lower conservative occlusion bound. When Re<Re 0 , occlusion is impossible, no matter how long the waves. In a real system, it is conceivable that wave coalescence events can produce local overshoots leading to occlusion even when Re<Re 0 . However, we have not observed this in our computations of spatially evolving films (see panels 16b and 17c). The value Re 0 being very low for all studied liquids, inertia is weak and significant overshoots are unlikely.

Based on the conservative occlusion bounds Re 0 and Re max , we delimit a region of conditional occlusion, Re 0 <Re<Re max (between the red dot-dashed and thick solid blue lines in panel 5b), where occlusion is theoretically possible, for long enough waves, but does not necessarily occur in a real system. That depends on whether the required waves actually develop in a spatially evolving film. And, in particular, whether the tube is long enough to accommodate the nonlinear wave dynamics required to produce them. In section 3.7, we will demonstrate that occlusion in this regime occurs through scenario II,

i.e. wave coarsening triggered by secondary instability, which shifts surface waves down the ascending branch of the dispersion curve toward increasingly low frequencies.

In panel 5a, the limit point LP is a turning point, where two solution branches meet.

As suggested by [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], we have checked their stability using transient periodic simulations with imposed Ψ started from selected travelling-wave solutions. In these computations, numerical noise is responsible for perturbing the flow. A full-fledged formal stability analysis, such as the one performed by [START_REF] Camassa | Traveling waves for a model of gravity-driven film flows in cylindrical domains[END_REF] based on the first-order model of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], is outside the scope of our current study.

We find that the lower solution branch is entirely unstable. As this branch does not connect to the flat-film primary flow d min =d max =1, the solutions thereon are thus highly unlikely in an experiment. Conversely, the upper branch is stable, at least w.r.t. periodic perturbations. However, wave trains of high-frequency travelling waves on falling liquid films are subject to secondary subharmonic and sideband instabilities that rely on wave interaction and lead to coalescence events [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF]. We have found that this triggers occlusion scenario II, which we will demonstrate in §3.7. Finally, when starting a periodic transient computation from a flat-film initial condition beyond the limit of travelling-wave solutions, the outcome is indeed occlusion, which we have checked (and validated with transient DNS) in appendix A (see panels 19a and 19b).

Panels 5c (case 2), 5d (case 3), and figure 6 (case 4) report the occlusion bounds Re max and Re 0 for the other three cases from table 1. In panel 5d, we have represented with open circles additional results obtained with our own travelling-wave DNS (see section 2.2), evidencing gratifying agreement with our WRIBL model predictions. Overall tendencies in these graphs are the same as for case 1. For case 4 (low-viscosity silicone oil), where the Kapitza instability is relevant, there is an additional feature. In panel 6a, the continuation curves of d min (Re,f =const) at high frequencies, e.g. the dot-dot-dashed curve at f /f max =1.2, exhibit two lobes that are each associated with a limit point. Wave profiles corresponding to these highlighted points LP1 and LP2 are represented in panel 6c. We see that the main wave hump is preceded by a small precursory ripple. Such socalled capillary ripples can produce multiplicity of solutions in falling liquid films, when they are close enough to interact with the following wave hump [START_REF] Kalliadasis | Falling Liquid Films[END_REF].

By contrast, the dashed curves in panel 6a, which pertain to low-frequency waves, each display only a single limit point (LP). Panels 6d and 6e show the corresponding wave profiles at those points. We see in panels 6c, 6d, and 6e that lower wave frequencies are associated with longer wavelengths, which we have implicitly assumed up until now.

Finally, panel 6f represents the wave profile corresponding to the limit point (LP) marked by a filled circle on the Re=2 curve in panel 6b. This curve lies in the region of conditional occlusion, Re 0 <Re<Re max (Re 0 =1.5, Re max =18.9). However, based on the wavelength of the waves in panel 6f, which span about 600 times the tube radius, it becomes clear that quite particular conditions can be required to cause occlusion in this regime. This underlines the conservative nature of the lower occlusion bound Re 0 .

In summary, due to the frequency dependence of travelling-wave solutions demonstrated in panels 5b, 5c, 5d, and 6b, three characteristic regimes can be delimited: (i)

impossible occlusion, when Re<Re 0 , (ii) conditional occlusion, when Re 0 <Re<Re max , and (iii) certain occlusion, when Re>Re max . In section 3.7, we will characterize these regimes with spatio-temporal computations and demonstrate the associated occlusion scenarios I and II. On the other hand, the frequency dependence of travelling-wave solutions opens the possibility of using coherent inlet forcing to prevent occlusion. By forcing regular waves of sufficiently high frequency, occlusion can be avoided versus a noise-driven wave evolution. We demonstrate this idea in section 3.8 based on spatiotemporal computations. In the next sections, 3.3 to 3.6, we first proceed to establish how different physical effects, i.e. gravity, axial viscous diffusion, inertia, and the core gas flow, affect the upper conservative occlusion bound Re max .

Role of gravity

The upper occlusion bound Re max can be significantly delayed by increasing the relative strength of gravity versus capillarity. We have quantified this by varying the Bond number Bo=ρ l g R ⋆2 /σ over several orders of magnitude at constant Laplace number La=6.3×10 -4 . This value of La corresponds to the high-viscosity silicone oil (case 1), where inertia is always negligible. Panel 7a represents travelling-wave solutions for different values of Bo. On each curve, we have imposed the frequency f max of the most amplified spatially growing linear waves.

All curves exhibit a limit value for the liquid Reynolds number Re. For Bo ≪ 1, where the Plateau-Rayleigh instability dominates, this limit is dictated by the onset of absolute instability (AI), as discussed in panels 2c and 2d. Beyond this point, prevailing waves are no longer characterized by the most-amplified spatially growing linear waves and f max is no longer defined. The solution curves stop abruptly. Occlusion in this regime is caused by absolute instability of the primary flow. This regime is not the focus of our study.

We focus on Bo 1, for which the solution curves in panel 7a are bounded by a limit point (LP), as discussed in panel 5a. This limit point marks the upper occlusion bound Re max , which increases significantly with increasing Bo, while the minimal core radius d min at LP significantly diminishes. The first observation is in line with the experiments of [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], who reported Re at the experimental occlusion onset for different liquids and tube radii but had no access to the wave dynamics. The authors also outlined an approach for predicting occlusion based on travelling waves, which relied on the conjecture that the wave height always reaches the tube radius as Re is increased toward the occlusion limit. Based on our results in panel 7a, this conjecture only holds when Bo is extremely large. Otherwise, travelling-wave solutions are lost well before their height attains the tube radius.

In panel 7b, we have replotted our results from panel 7a in terms of the liquid hold-up V l /π/Λ. We see that the maximally achievable liquid hold-up increases significantly with Bo. This again is in line with the experiments of [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], who had observed that the mean film thickness h at occlusion is significantly greater than in the gravity-free Plateau-Rayleigh configuration [START_REF] Gauglitz | An extended evolution equation for liquid film breakup in cylindrical capillaries[END_REF].

The variation of Bo in figure 7 amounts to varying the gravitational acceleration for constant liquid and tube radius. Of course, it is either very difficult (low Bo range, microgravity experiments) or impossible (high Bo range) to realize the displayed range of variation 0.1 Bo 100 practically. A more practicable way to quantify the desired effect is to increase the tube radius R for a constant liquid. This amounts to increasing Bo at constant Ka=σ ρ -1 l g -1/3 µ -4/3 l . We have performed such computations for Ka=3.3×10 -3 , corresponding to the high-viscosity silicone oil (case 1). Results are represented in panels 8a and 8b, and these display the same overall behaviour as figure 7.

The delay of Re max results from a gravity-induced streamwise distortion of surface waves, which favours variations in (stabilizing) axial interface curvature over variations in Bo. This allows to counter a larger variation in (destabilizing) azimuthal curvature, and, thus, to saturate waves of greater amplitude.

The saturation mechanism was identified by [START_REF] Frenkel | Annular flows can keep unstable films from breakup: nonlinear saturation of capillary instability[END_REF] for pressure-driven core annular flows, where it is due to a viscous interfacial shear stress, and demonstrated experimentally for falling liquid films on cylindrical fibres by [START_REF] Quéré | Thin films flowing on vertical fibers[END_REF]. In both studies, the saturation took hold in the weakly-nonlinear regime, limiting waves to a very small amplitude. Thus, it was perceived as suppressing macroscopically-visible waves altogether. In our case, the mechanism takes hold in the strongly-nonlinear regime, allowing to saturate travelling waves of greater amplitude as Bo is increased.

The role of gravity is kinematic, i.e. it makes initial surface elevations travel faster than surface depressions [START_REF] Dietze | On the Kapitza instability and the generation of capillary waves[END_REF], and this produces the distortion of wave fronts observed in panel 8c. We point out that this effect was precluded in [START_REF] Jensen | Draining collars and lenses in liquid-lined vertical tubes[END_REF], where surface waves were constructed in the form of perfectly symmetrical unduloids.

Role of axial viscous diffusion

Axial viscous diffusion, emanating from the ∂ xx u term in (2.1), is represented in our model (2.5) through the terms with coefficients J j , K j , L j , and M j . These terms appear when the underlying long-wave expansion is performed up to second order. This was not done in the previous modelling works of [START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF], [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF][START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 2. air-driven flow[END_REF], [START_REF] Zhou | Viscoelastic liquid film flowing down a flexible tube[END_REF], [START_REF] Liu | Stability of viscous film flow coating the interior of a vertical tube with a porous wall[END_REF], and [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF]. So, the question is whether retaining these terms is worth the effort. We have checked this in panels 9a, 9b, and 9d, representing travelling-wave solutions for cases 1, 3, and 4. The three cases span a wide range of the capillary number Ca=µ l U/σ, which relates viscous to capillary stresses, and which quantifies the relevance of axial viscous diffusion.

In figure 9, solid lines correspond to the full model (2.5) and dot-dashed lines to the limit J i =K i =L i =M i =0. Further, we distinguish two types of continuations. On the blue curves, which contain no symbols, we have fixed the wavelength Λ=12.56. On the black curves with symbols, we have imposed the spatially most amplified frequency of linear waves f =f max . In that case, axial viscous diffusion may affect the solution both nonlinearly and linearly (via wave selection), but the latter effect is negligible for the cases considered here. By comparing the solid and dashed curves in panels 9a, 9b, and 9d, we see that axial viscous diffusion precipitates occlusion in general and Re max in particular and this effect increases with Ca.

For the high-viscosity silicone oil (panel 9a, case 1: Ca=0.37), accounting for axial vis- For case 3 (panel 9b, Ca=0.028), which corresponds to experimental run 13 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], axial viscous diffusion dictates even the nature of the occlusion mechanism. For the continuation at f =f max (black lines with symbols in panel 9b), occlusion is caused by absolute linear instability (AI) when axial viscous diffusion is neglected (dot-dashed line), and by nonlinear loss of travelling-wave solutions (LP) when it is taken into account (solid line).

For the lowest Ca (panel 9d, case 4: Ca=0.012), neglecting axial viscous diffusion does not meaningfully affect Re max . It is however responsible for the formation of a second lobe on the black dot-dashed curve marked by an asterisk, which results from a change in the number of precursory ripples, as discussed w.r.t. panels 6a and 6c.

The nonlinear mechanism by which axial viscous diffusion precipitates occlusion is

Ψ xx = Fr 2 Re ∂ x F xx π(1 -d 2 )
, F xx = 2π Such waves are known as roll waves and the grey area in panel 9c delimits their existence.

The limit point on the solid black line (filled circle) lies within this region, while the limit point on the dot-dashed line (asterisk) does not. This is illustrated by the corresponding streamline patterns in figure 11. Also, we see that the roll wave in panel 11a conveys with it a large gas bubble containing a toroidal vortex. This particular flow pattern results from the aerostatic pressure difference imposed here, i.e. Ψ=1. Other gas flow scenarios will be discussed in section 3.6. Finally, we point out that the lower branches in panel 9 are the stable ones, and that only the full model, which accounts for axial viscous

diffusion, predicts roll wave solutions on these branches.

Role of inertia

Inertia constitutes the growth mechanism of the Kapitza instability, which is responsible for the formation of surface waves on planar falling liquid films, and which becomes relevant in our current problem when Re is sufficiently large. This is the case for the lowviscosity silicone oil film (case 4). In figure 12, we have represented nonlinear travellingwave solutions (panel 12a) and dispersion curves of the linear spatial growth rate α i (panel 12b) for this case. Solid curves were obtained with our full model (2.5), and dashed curves by neglecting inertia, i.e. by setting S i =F ij =G ij =0 in (2.5) and (2.11).

In panel 12a, we confront two types of travelling-wave continuations. On the blue lines without symbols, we have fixed the wavelength Λ=12.56. In this case, the effect of inertia is purely nonlinear. We see that this effect significantly advances the limit point of the solid blue versus the dashed blue curve (by approximately 20%). Inertia increases the cumulated strength of the destabilizing mechanisms (Plateau-Rayleigh and Kapitza instabilities) versus the stabilizing ones (capillarity and viscous dissipation). This limits the ability to saturate nonlinear travelling waves to smaller amplitudes, i.e. lower Re.

On the black lines with symbols, we have imposed f =f max . Here, the effect of inertia Re max . When it is neglected, occlusion occurs to to absolute instability (AI).

In addition to shifting Re max , inertia also changes the nature of the occlusion bound in panel 12a. If inertia is neglected, Re max is dictated by absolute instability (AI, asterisk), otherwise, it is dictated by a nonlinear turning point (LP, filled circle). This qualitative change was also observed for the role of axial viscous diffusion in figure 9.

In the case of the classical Plateau-Rayleigh configuration, Bo → 0, inertia only affects the dynamics of the film's evolution [START_REF] Dietze | Films in narrow tubes[END_REF] but not the occlusion limit itself, because the latter is determined by a static equilibrium state [START_REF] Everett | Model studies of capillary condensation[END_REF].

Role of core gas flow

The gas in the core of the cylindrical tube can affect the liquid film falling down its inner wall in two ways: (i) through the normal and tangential interfacial viscous stresses τ n g and τ t g , which are scaled by Π µ in (2.3), and (ii) through the interfacial gas pressure p g , which is scaled by Π ρ in (2.3). Up until now, we have assumed that the gas is subject to an aerostatic pressure difference, Ψ=1, meaning that gravity is fully compensated in the gas. In that case, the gas flow is enslaved to the hydrodynamics of the liquid film and both the above effects are negligible. Nonetheless, this one-way coupling produces intricate flow structures in the gas, as has been suggested by panel 11a and will be discussed further based on panel 14a.

Conversely, when the gas flows counter-currently at a Reynolds number Re g of sufficient magnitude, it does affect the hydrodynamics of the liquid film, and this two-way coupling can significantly precipitate the upper occlusion bound Re max . We show this in figure 13 based on case 4, which produces surface waves of sufficient amplitude to significantly constrict the core gas flow. Panel 13a represents travelling-wave solutions obtained by continuing Re at f =f max for different gas flow configurations. The solid black curve with a filled circle at the limit point corresponds to Ψ=1, and the dashed curve to a countercurrent gas flow with Re g =-17.27. Comparing these two curves, we see that Re max is significantly precipitated due to the counter-current gas flow (by roughly 25%). We have checked that the gas flow does not meaningfully affect the most amplified frequency f max , and thus this precipitation is a nonlinear effect. The dashed curve also displays two qualitative changes. First, its lower branch turns forward at a second, lower, limit point, producing travelling-wave solutions of constant amplitude that exist far beyond Re max . However, as we show in appendix C, these solutions are inaccessible in a real system. Thus, Re max remains a representative occlusion bound also in the counter-current configuration. Second, solutions on the portion of the curve between the asterisk and the new limit point are stable, in contrast to those on the lower branch of the solid black curve, where Ψ=1.

Panel 13b illustrates the effect of the gas flow in another way. Here, Re g is varied at a fixed liquid volume V l , while maintaining f =f max . Such a scenario may be relevant for mucus films in pulmonary capillaries. We see that travelling-wave solutions are lost when the magnitude of the counter-current gas flow |Re g | is increased beyond a limit point.

The effect of the gas flow on the occlusion bound Re max observed in panel 13a results from the pressure coupling between gas and liquid and not from the gaseous interfacial viscous stresses. This is evidenced by the curves without symbols in panel 13a, which pertain to different limits of the full-model continuation at Re g =-17.27 (dashed curve).

The blue solid curve without symbol pertains to the limit Π µ =0 and the dot-dashed curve to the limit Π ρ =0. Deactivating the gaseous viscous stresses (Π µ =0, solid blue curve without symbols), changes Re max very little. By contrast, deactivating the gas pressure effect (Π ρ =0, dot-dashed curve), almost entirely negates the gas-induced precipitation of Re max versus the aerostatic reference case (solid black curve with filled circle).

To elucidate the role of the gas pressure p g , panels 13c and 13d represent wave profiles and corresponding profiles of the pressure gradient ∂ x p g for the two marked limit points (LP) in panel 13a. Thereby, ∂ x p g is determined from (2.5c) and corresponds to the gas pressure gradient at the film surface ∂ x p g | r=d . We have dropped the accompanying subscripts for convenience. In panel 13d, ∂ x p g is multiplied with Fr 2 , in order to normalize the actual pressure gradient with the aerostatic one. Grey/white zones between the dashed lines in panels 13c and 13d distinguish regions of negative/positive ∂ x p g in the gas. In addition, figure 14 represents streamlines in the wave-fixed reference frame corresponding to the wave profiles in panel 13c.

Focussing on the counter-current case, panel 14b, we see that the wave hump produces a constriction through which the gas (flowing from right to left) must pass, similar to a G. F. Dietze, G. Lavalle and C. Ruyer-Quil de Laval nozzle that first accelerates and then decelerates the flow. This produces strong positive/negative pressure gradients to the right/left of the wave maximum (see dashed line in panel 13d) and these tend to push the liquid within the film toward the wave hump (arrows in panel 13c). This mechanism tends to increase the wave amplitude for a

given Re and thus to advance the occlusion bound Re max . It requires a sufficiently strong relative motion between the traveling wave and the gas. This does not necessarily require an upward gas flow. Indeed, the diamond-marked curve in panel 13b shows that occlusion may be brought about even when the gas is co-current, at least if |Re g | is sufficiently small.

In the aerostatic case, Ψ=1, axial pressure variations in the core are very weak (solid line in panel 13d) and the gas has no incidence on the occlusion limit. Vice versa, the flow structure in the gas is strongly affected by the liquid film. We see this in panel 14a, which pertains to case 4 (low-viscosity silicone oil), and also in panel 11a, which pertains to case 1 (high-viscosity silicone oil). In both examples, the gas flow is enslaved to the liquid surface flow (through equation 2.3a), which exhibits interfacial stagnation points around the wave hump, as a result of the moving-frame vortex in the liquid film. Due to this, the gas flow is divided into an intricate system of vortices. The constellation of these vortices is dictated by the strength of the liquid vortex in the wave hump.

Occlusion scenarios in a spatially evolving film

In section 3.2, we introduced the conservative occlusion bounds Re 0 and Re max based on the limit of travelling-wave solutions and delimited three regimes for a naturally evolving wavy falling liquid film: (i) impossible occlusion, when Re<Re 0 , (ii) conditional occlusion, when Re 0 <Re<Re max , and (iii) certain occlusion, when Re>Re max . We now demonstrate these regimes based on spatio-temporal computations with our model (2.5),

where we have mimicked naturally-evolving waves through the noisy inlet condition

(2.6b), using ǫ 2 =10 -5 . For the gas, we have imposed Ψ=1.

We start by reproducing two of the occlusion experiments of [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF],

i.e. their experimental runs 20 and 13 which identified the experimental occlusion point at Re=0.0497 and Re=0.258, respectively. These two runs correspond to our cases 2 and 3 in table 1, for which we have determined the lower/upper occlusion bounds in §3.2 as Re 0 =0.012/Re max =0.038 (panel 5c), and Re 0 =0.075/Re max =0.39 (panel 5d), respectively. Based on this, experimental run 20 (our case 2) lies in the region of certain occlusion and experimental run 13 (our case 3) in the region of conditional occlusion. This is confirmed by panels 15a and 15b, which represent travelling-wave solutions continued in terms of the wave frequency f (right abscissae) and dispersion curves of the linear growth rate α i (left abscissae) for these two runs.

The experiments of [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] were performed in a long vertical tube equipped with a conductance probe that allowed to detect occlusion locally. However, no information on the film's surface dynamics upstream of this point was available.

We provide this missing information in panels 15c and 15d, which represent snapshots of our spatio-temporal computations reproducing experimental runs 20 and 13. These computations were performed on an L ⋆ =1 m domain. We see that the two cases reach occlusion by two very different routes, as a result of distinct wave dynamics.

For experimental run 20 (our case 2), most of the dispersion curve in panel 15a (grey area) lies beyond the limit point of travelling wave solutions (f max <f LP ). Thus, the most amplified surface waves emerging from the linear regime don't possess any saturated travelling state and inevitably cause occlusion by extending their growth into the nonlinear regime (panel 15c). This is occlusion scenario I, which we had already observed for case 1 (high-viscosity silicone oil) in figure 4. Supplementary movie M2 shows the computation from panel 15c in action. We remind the reader that the unphysical plug reopening observed therein (and in supplementary movie M3) is due to the crude numerical core radius limitation introduced in section 2.1. First results in appendix E suggest that this artefact is negated by our improved WRIBL model (figure 22).

By contrast, for experimental run 13 (our case 3), the most amplified linear waves lie well within the range of travelling-wave solutions (f max >f LP in panel 15b). Occlusion must thus occur through a different route, as evidenced by the wave profile in panel 15d. In the first part of the tube (x 130), a quite regular train of saturated-amplitude travelling waves develops. These result from linear wave selection and their minimal core radius d min is far from reaching the tube axis. However, individual waves within the wave train are quite narrowly spaced and thus prone to the secondary subharmonic/sideband instabilities typically associated with single-humped waves in planar falling liquid films [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF][START_REF] Chang | Nonlinear evolution of waves on a vertically falling film[END_REF]. These instabilities result from wave interactions triggered by residual noise and lead to wave coalescence, which is the first step in the so-called wave coarsening dynamics identified by Chang et al. (1996b). This cascaded dynamics produces ever longer, more dangerous, waves by successive coalescence/absorption events and thus represents an inherent route toward occlusion. The arrow in panel 15d highlights one such coalescence event, which produces a wave of greater amplitude that will occlude the tube and subsequently absorb several smaller waves downstream. This is occlusion scenario II, which is intermittent because the initiating wave coalescence is occasional, in contrast to scenario I (panel 15c), where (almost) every wave grows to occlude the tube. Due to this intermittence, gas bubbles resulting from successive occlusion events are much more unevenly distributed. Supplementary movie M3 shows the computation from panel 15d in action.

We now focus on case 3 and compute the spatio-temporal wave evolution for values of Re below the experimental occlusion threshold Re=0.258 determined in run 13 of [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF]. This experimental threshold lies well within the region of conditional occlusion, which is bounded by Re 0 =0.075 and Re max =0.393 (panel 5d).

Panel 16a represents a spatio-temporal computation at a lower Reynolds number Re=0.14>Re 0 , which still lies in the conditional occlusion regime. This computation was performed on a domain of L ⋆ =2 m in length, i.e. twice longer than the one in panel 15d. We see that occlusion still occurs, through scenario II, but that the occlusion point has shifted greatly downstream w.r.t. panel 15d. The reason for this is that the frequency/wavelength of the limiting travelling waves greatly decreases/increases with decreasing Re (see panel 5d) and, thus, much more space is needed for these to develop from the coarsening dynamics. The downstream shift in the occlusion point explains why [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], who checked for occlusion in a fixed region, found an experimental threshold for Re that is greater than the lower conservative bound Re 0 . It also clearly shows that occlusion in the conditional regime (through scenario II) depends on whether the tube is sufficiently long to accommodate formation of those waves that are able to cause occlusion. As these waves may be much longer than those emerging from linear selection, occlusion may require unrealistically long tube lengths, depending on Re. Supplementary movie M4 shows the computation from panel 16a in action.

Conversely, occlusion in an infinitely long tube will eventually always occur if Re lies within the conditional regime, at least for a naturally evolving film. This is due to the nature of the coarsening dynamics, which halts only after truly solitary waves have formed. Thus, to avoid occlusion in an infinitely long tube, Re must be decreased below the lower conservative bound Re 0 , where even solitary waves are too small to cause Although occlusion is indeed avoided, surface waves are extremely small in amplitude, which thwarts their beneficial effect on heat and mass transfer [START_REF] Yoshimura | Enhancement of mass transfer into a falling laminar liquid film by two-dimensional surface waves-some experimental observations and modeling[END_REF][START_REF] Albert | Direct numerical simulation of interfacial mass transfer into falling films[END_REF]. This underlines the conservative nature of the lower occlusion bound Re 0 . We will show in the next section that there is another way to prevent occlusion, namely through coherent inlet forcing in the conditional occlusion regime, which allows maintaining surface waves of significant amplitude.

Finally, figure 17 shows snapshots of three representative spatio-temporal computations for case 4 (low-viscosity silicone oil), where inertia, which drives the the Kapitza instability, is relevant. Panels 17a, 17b, and 17c correspond to the regimes of certain, conditional, and impossible occlusion. These display the same main features observed for the high-viscosity liquids (panels 15c, 16a, and 16b), only that the length scales of the occlusion processes are much shorter.

Coherent inlet forcing to prevent occlusion

For Re>Re 0 , travelling-wave solutions are bounded by a limiting frequency f LP , below which they cannot exist and occlusion is inevitable (see panels 5d and 6b). Conversely, by increasing f above f LP , travelling-wave solutions are recovered and occlusion can be avoided. We thus test the idea of forcing high-frequency waves through coherent inlet forcing in a spatially-evolving falling liquid film, in order to prevent occlusion without having to reduce Re below the conservative bound Re 0 , where surface waves all but disappear (see panels 16b and 17c).

This idea works only in the conditional occlusion regime, Re 0 <Re<Re max . Indeed, in the regime of certain occlusion, Re>Re max , where f LP >f max , one would have to force waves of frequency f >f max . Such short waves are very sensitive to the secondary instabilities causing wave coalescence in closely-packed wave trains [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF], and rapidly loose the signature of the inlet forcing. By contrast, in the conditional occlusion regime, f LP can be low enough (the lower Re, the lower f LP ) to allow forcing a wave train that is not subject to secondary instability.

We demonstrate this by running the spatio-temporal computations in panels 16a (case 3) and 17b (case 4) again with additional coherent inlet forcing. For this, we activate the monochromatic perturbation in (2.6b) with a relative amplitude ǫ 1 =0.1. Meanwhile, the inlet noise remains active and its level unchanged ǫ 2 =10 -5 . We stress that ǫ 2 in (2.6b) is the scale factor of a Fourier series which 1000 terms. Thus, the actual noise amplitude is much greater than suggested by the value of ǫ 2 . It is about 10% of ǫ 1 for the runs presented here. The forcing frequency f is chosen greater than the limit value f LP of travelling-wave solutions, as obtained from panels 5d (case 3, Re=0.14, f LP =0.41f max ) and 6b (case 4, Re=5, f LP =0.35f max ), respectively. Panels 18a and 18b show the results of our spatio-temporal computations with additional coherent inlet forcing and are to be compared directly to panels 16a and 17b, where the wave dynamics is purely noise-driven. We see that the additional coherent inlet forcing produces a regular train of saturated-amplitude travelling waves that do not occlude the tube. This wave train is maintained over the entire length of the domain and thus robust w.r.t. to secondary instability. In contrast to panels 16b and 17c, where occlusion was avoided by reducing Re below the lower conservative bound Re 0 , the falling films in panels 18a and 18b display quite substantial surface waves, expected to be much more beneficial for inter-phase heat and mass transfer.

Conclusion

We have studied numerically the occlusion of a narrow vertical cylindrical tube by an axisymmetric wavy falling liquid film in contact with a laminar core gas flow. We have focused on conditions where the effect of gravity is at least comparable to that of capillarity, i.e. Bo 1, and the gas Reynolds number Re g is quite low. In that limit, occlusion is caused by spatially growing waves and not by absolute linear instability. [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF], while scenario II (panel 15d, movie M3) applies to run 13 in [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF].

We have delimited scenarios I and II based on an upper conservative occlusion bound Re max , which corresponds to the limit point of travelling-wave solutions at the spatially most amplified frequency f max of linear waves. And, a lower conservative bound Re 0 , below which travelling-wave solutions always exist, no matter how great the wavelength [START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF]. Determining Re max requires computing f max and thus we have simultaneously solved the linear stability problem in our numerical continuations. This had not been attempted in previous studies [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF][START_REF] Ding | Thermocapillary effect on the dynamics of liquid films coating the interior surface of a tube[END_REF].

Based on the upper and lower conservative bounds Re max and Re 0 , as well as spatiotemporal computations, we have delimited three possible regimes for a naturally evolving wavy film: (i) certain occlusion: Re>Re max . The most amplified surface waves, which typically emerge in an experiment, do not possess travelling states. Occlusion is inevitable and occurs through scenario I; (ii) conditional occlusion: Re 0 <Re<Re max . Occlusion is theoretically possible through scenario II (for long enough waves), but does not necessarily occur in a real system. That depends on whether the coarsening dynamics can produce sufficiently long waves and whether the tube is long enough to accommodate this (panels 15d and 16a, movies M3 and M4); (iii) impossible occlusion: Re<Re 0 . All possible surface waves, no matter how long, are safe.

We have shown that occlusion can be actively prevented in the regime of conditional occlusion (Re 0 <Re<Re max ) by forcing waves of sufficiently high frequency through coherent inlet forcing (panels 18a and 18b). This allows to maintain surface waves of nonnegligible amplitude, in contrast to the regime of impossible occlusion (Re<Re 0 ) where waves are almost invisible (panels 16b and 17c). Given the well documented positive effect of surface waves on inter-phase heat/mass transfer [START_REF] Yoshimura | Enhancement of mass transfer into a falling laminar liquid film by two-dimensional surface waves-some experimental observations and modeling[END_REF], coherent inlet forcing is an attractive route to enhance transfer while avoiding occlusion.

We have computed the lower and upper conservative bounds Re 0 and Re max for the four working liquids considered here (figures 5 and 6). Further, we have quantified how key parameters (Bo and Re g ) as well as several physical effects (axial viscous diffusion and inertia) affect the upper conservative bound Re max , which delimits the regimes of certain and conditional occlusion. We have found that gravity greatly delays this bound, significantly increasing Re max as the Bond number Bo is increased (see panel 8a). Unless

Bo is very large, the height of travelling waves at Re=Re max is far from reaching the tube radius. In contrast to what [START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF] conjectured, such solutions are thus lost abruptly and not due to the wave height reaching the tube radius continuously.

We have found that a laminar counter-current gas flow can significantly precipitate occlusion for low-viscosity liquids (low-viscosity silicone oil, case 4), where the effect of the gas pressure becomes relevant in the liquid force balance. At the strongest gas flow studied (Re g =-17), we have observed a 25% reduction of Re max (panel 13a) versus a situation where the gas is quiescent. We have found that the gas affects the liquid film mainly through inter-phase pressure coupling. Surface waves locally constrict the gas flow, which is accelerated and then decelerated as it passes through (panel 14b). This produces a pressure minimum above the wave crest that sucks liquid toward the wave, tending to increase its amplitude and, thus, the danger of occlusion.

For high-viscosity liquids (high-viscosity silicone oil and aqueous glycerol solution, cases 1 and 3), axial viscous diffusion greatly precipitates occlusion, i.e. greatly reduces Re max (by 74% in panel 9a). This may be relevant for mucus films within the first five generations of the human respiratory network, where the capillary radius is large enough for gravity to be relevant [START_REF] King | Rheological properties of microliter quantities of normal mucus[END_REF][START_REF] Lewis | Quantification of airway diameters and 3d airway tree rendering from dynamic hyperpolarized 3He magnetic resonance imaging[END_REF][START_REF] Grotberg | Respiratory fluid mechanics[END_REF]).

For low-viscosity liquids (low-viscosity silicone oil, case 4), we have observed that Re max becomes large enough for inertia to be relevant (Re max =18.9 in panel 6a). This activates the Kapitza instability, which tends to increase the amplitude of travellingwaves, increasing the danger of occlusion. We have found that inertia significantly precipitates occlusion, decreasing Re max versus an inertialess computation (by 20% in panel 12a). This may be relevant for falling film micro-reactors [START_REF] Seebauer | Tube bundle falling film microreactor for performing gas liquid reactions[END_REF].

Finally, we have found that neglecting axial viscous diffusion (high-viscosity liquids, panel 9a) or inertia (low-viscosity liquids, panel 12a) can change the nature of the occlusion mechanism, from being dictated by the loss of travelling-wave solutions to being dictated by absolute linear instability.

Our transient computations on long open domains have been performed with a crude numerical core radius limitation allowing to continue these past occlusion events. Although this technique does not affect the wave dynamics leading up to occlusion, which is the focus of our current study, it does not represent the physics of liquid plugs to satisfaction. To remedy this, we have introduced an improved version of our WRIBL model (2.5) by augmenting it with an additional force term (4.1), allowing to form stable travelling pseudo-plugs (figure 22). First results reported in appendix E suggest that this improved model negates the main artefact of the numerical core radius limitation, i.e. the reopening of liquid plugs. Confronting our improved model with full-fledged plug models [START_REF] Ubal | Stability of the steadystate displacement of a liquid plug driven by a constant pressure difference along a prewetted capillary tube[END_REF][START_REF] Suresh | The effect of gravity on liquid plug propagation in a two-dimensional channel[END_REF]) is an enticing prospect for future work. for the unknowns U=(h, D ξ h, D ξξ h, a i , D ξ a i , D ξξ a i , b i-1 ) ∀ 1 i n. Inverting (4.5) leads to an autonomous dynamical system of dimension 4n + 3. This dynamical system was solved with the continuation software AUTO07P [START_REF] Doedel | Continuation and bifurcation software for ordinary differential equations[END_REF], using a predictorcorrector method [START_REF] Kalliadasis | Falling Liquid Films[END_REF]. Our computations were performed with η=10 -6 and n = 20, which we have determined based on convergence tests.

Appendix C: Travelling waves on the lower branch in panel 13a

We discuss further the travelling-wave solutions on the lowest branch of the dashed curve in panel 13a, which extends beyond the limit point Re=Re max , as the result of interfacial viscous stresses exerted by the gas. Panel 20a replots that curve (dashed line)

and the corresponding curve for the aerostatic configuration (solid line) in terms of the pressure gradient Ψ=∆p g Fr 2 . The filled circles therein correspond to the limit points (LP) in panel 13a. The segment of the dashed curve that we are interested in lies above the asterisk, and we have checked the stability of solutions thereon, using transient periodic computations with imposed Ψ (we recall that the red portion of the curve, between circle and asterisk, is stable). As Re is increased on this curve segment, solutions go from being stable to being subject to an oscillatory instability, similar to that found in planar channels (Lavalle et al. to be published). In both situations, our transient periodic computations do not produce occlusion.

However, these states are inaccessible in a real system, where Re is fixed and the counter-current gas flow is imposed through Ψ. The arrow in panel 20a indicates how Ψ would need to be varied to attain the desired states, starting from the aerostatic situation (marked by a cross). Besides that this would require a two hundred fold increase of Ψ, it turns out that there is no continuation path (at Re=const) from the solid curve, where -17.27)=0.73, while gradually increasing Ψ from Ψ=1 according to a sigmoid function.

Ψ=1
This computation produces occlusion long before Ψ has reached the target value Ψ=237, which corresponds to the upper solution branch in panel 20a. We may thus conclude that Re max , which is based on the limit point (LP) in panel 13a, remains a representative occlusion bound also in the counter-current configuration.

Appendix D: Effect of axial viscous diffusion in panel 4a

Figure 21 represents a version of the computation from panel 4a where we have deactivated axial viscous diffusion, by setting J j =K j =L j =M j =0 in (2.5b). In contrast to panel 4a, no occlusion is observed over the entire length of the tube, which is in qualitative contradiction with the experiments of [START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF].

Appendix E: Improved representation of liquid plugs

Inspired by the use of an additional attractive-repulsive force term to model contact line problems within the framework of film models [START_REF] Thiele | Sliding drops in the diffuse interface model coupled to hydrodynamics[END_REF], we add a repulsive azimuthal capillary term to the RHS of our model equation (2.5b):

Π θ = We -1 exp λ 1 - d(x, t) d plug ∂ x κ θ , κ θ = - 1 d 1 - 1 2 (∂ x d) 2 . (4.1)
As a result of its exponential variation with d (which can be scaled with λ), the effect of this term is significant only in the vicinity of the predefined radius d plug , which is set to a small fraction of the tube radius R. At d=d plug , Π θ exactly compensates the capillary pressure gradient due to (destabilizing) azimuthal curvature. Furthermore, by adequately choosing d plug and λ, the primary flow with a cylindrical film surface at d=d plug can be rendered linearly stable at all wave numbers (we have checked this with spatial linear stability analysis). This allows representing liquid plugs as stable travelling liquid annuli, or pseudo-plugs, that almost entirely fill the tube cross section. , 30, 40, 50, 60, 70, 80, 90 (t ⋆ ν =ν 1/3 l g -2/3 =0.0037 s), during the formation of a pseudo-plug. The four last profiles have been shifted by one wavelength to better guide the eye; (b) fully developed pseudo-plugs: t ⋆ /t ⋆ ν =300. The periodic surface profile has been reproduced three times here.

be alleviated by implementing a locally refined grid in the plug regions, following the example of [START_REF] Lister | Capillary drainage of an annular film: the dynamics of collars and lobes[END_REF] for the representation of locally thin liquid films. We leave this task, and the confrontation of our pseudo-plug model with full-fledged plug models for pressure driven core-annular flows [START_REF] Ubal | Stability of the steadystate displacement of a liquid plug driven by a constant pressure difference along a prewetted capillary tube[END_REF] or inclined planar channel flows [START_REF] Suresh | The effect of gravity on liquid plug propagation in a two-dimensional channel[END_REF], to future work. In the latter case, our additional force term (4.1) would have to be adapted to a planar geometry.

Figure 1 .

 1 Figure1. Falling liquid film (subscript l) lining the inner surface of a narrow cylindrical tube in contact with a laminar gas flow in the core (subscript g). (a) Problem configuration and notations. The tube radius R ⋆ is the length scale, the star superscript designating dimensional quantities; (b) example of a travelling-wave computation with our model (2.5): case 4 in table 1 (low-viscosity silicone oil), Re=14.6, Reg=-17.27 (limit point marked by asterisk in panel 13a). Streamlines in the wall-fixed reference frame within the liquid (blue lines) and gas (red lines).

  1,2,3. The code employed for this is based on a Crank-Nicolson time discretization and first-order central differences for spatial discretization (Dietze & Ruyer-Quil 2015); (ii) continuation of periodic travelling wave solutions that are stationary in the wave-fixed coordinate ξ=x-c t moving with the wave celerity c. These computations were performed with the continuation software AUTO07P (Doedel 2008); (iii) transient open-domain computations with inlet and outlet

Figure 2 .

 2 Figure 2. Linear instability predictions for three cases from table 1: case 1 (dot-dashed lines, diamonds); case 3 (solid lines, open circles); case 4 (dashed lines, filled circles). The core fluid is subject to an aerostatic pressure gradient Ψ=1. Lines correspond to the model dispersion relation (2.11). (a) Temporal growth: α∈ R, ω=ωr+i ωi in (2.10). Symbols correspond to Orr-Sommerfeld solution, following the work of Hickox (1971); (b) spatial growth: ω ∈R, α=αr+i αi in (2.10). Crosses highlight maximal spatial growth; (c) dispersion curves under increasing Re∈ [0.1, 0.54]. Blue line traces maximum growth rate αi(ωmax); (d) most amplified angular velocity ωmax versus Re. Asterisks mark absolute linear instability threshold (AI).

Figure 4 .

 4 Figure 4. Nature of the prevailing surface waves that cause occlusion. Spatio-temporal computation reproducing the second version of the experiment in panel 3a of Camassa et al. (2014), which was performed in a longer L ⋆ =1 m tube: case 1 in table 1, Re=2.3×10 -4 . (a) Snapshot of the computed film profile, showing individual surface waves that cause occlusion in the region x 130; (b-d) frequency spectra at different streamwise locations leading up to the occlusion region (marked by dashed lines in panel a): (b) x=20; (c) x=70; (d) x=110.Vertical blue lines at f /fmax=1.12 mark limit of travelling wave solutions determined by numerical continuation (see panel 6b). fmax designates spatially most amplified frequency of linear waves.

  amplitude. Consequently, low-frequency travelling waves are more prone to occlusion than high-frequency waves. This effect was demonstrated by Camassa et al. (2016) and Ding et al. (2019) with a first-order asymptotic model and a lubrication model, respectively.

Figure 5 .

 5 Figure5. Limits of travelling-wave solutions obtained from numerical continuation with our model (2.5). Cases 1 (panels a, b), 2 (panel c), and 3 (panel d ) in table 1. (a) Continuation of Re at wave frequency f /fmax=0.5 (long dashes), 0.8 (short dashes), 1 (solid), 1.1 (dot-dashed), and 1.2 (dot-dot-dashed). Vertical blue line marks limit point (LP) for spatially most amplified linear waves, where f =fmax and Re=Remax; (b-d) frequency dependence at Re=const. Dashed blue: linear spatial growth rate αi at Re=Remax (right abscissae), other: travelling-wave solutions at Re=const (left abscissae). (b) Case 1: Re=Re0=6 × 10 -5 (dot-dashed red), Re=7.5 × 10 -5 , 1 × 10 -4 (thin solid), Re=Remax=1.91 × 10 -4 (thick solid blue), and Re=2.3 × 10 -4 (dotted with asterisk); (c) case 2: Re=Re0=0.012 (dot-dashed red), Re=0.016, 0.025 (thin solid), and Re=Remax=0.038 (thick solid); (d) case 3: Re=Re0=0.075 (dot-dashed red), Re=0.1, 0.14, 0.2 (thin solid), and Re=Remax=0.393 (thick solid). Open circles correspond to data obtained with our own travelling-wave DNS (see section 2.2).

Figure 6 .

 6 Figure 6. Limits of travelling-wave solutions for case 4 in table 1. (a) Continuation of Re at constant wave frequency f /fmax=0.5 (long dashes), 0.8 (short dashes), 1 (solid), 1.1 (dot-dashed), and 1.2 (dot-dot-dashed); (b) frequency dependence at Re=const. Dashed: linear spatial growth rate αi at Re=Remax=18.9 (right abscissa); other: travelling-wave solutions (left abscissa) at Re=Re0=1.5 (dot-dashed red), Re=2, 5, 10 (thin solid), and Re=Remax (thick blue); (c-e) film profiles for LPs marked by filled circles in panel a; (c) f =1.2fmax. Solid: LP1; dashed: LP2; (d ) f =0.8fmax; (e) f =0.5fmax; (f ) profile for marked LP in panel b.

Figure 7 .

 7 Figure7. Delay of occlusion by increasing the Bond number Bo=ρ l g R ⋆2 /σ at constant Laplace number: La=σ ρ l R ⋆ /µ 2 l =6.27×10 -4 (case 1). Travelling-wave solutions at f =fmax. Symbols mark onset of absolute linear instability (AI) or limit points (LP). The limit points correspond to the upper occlusion bound Remax. Circles: Bo=0.1; squares: Bo=1; asterisks: Bo=11.1 (case 1); crosses: Bo=50; diamonds: Bo=100. (a) Minimal core radius dmin versus liquid Reynolds number; (b) minimal core radius versus liquid holdup.

Figure 8 .

 8 Figure 8. Delay of occlusion by increasing the Bond number Bo at constant Kapitza number: Ka=σ ρ -1 l g -1/3 µ -4/3 l =3.29×10 -4 (case 1). This amounts to increasing the tube radius R ⋆ while maintaining constant fluid properties. See figure 7 for symbol attributions. (a) Minimal core radius dmin versus liquid Reynolds number; (b) minimal core radius versus liquid holdup; (c) surface profiles at AI and LP; (d) corresponding profiles of axial curvature ∂xxd.

(

  destabilizing) azimuthal curvature. Panel 8c illustrates this distortion, which is marked by a compression of the leading wave front and an elongation of the trailing front. The different surface profiles represented correspond to the absolute instability thresholds (AI) and limit points (LP) highlighted in panel 8a. Panel 8d represents corresponding profiles of the axial curvature ∂ xx d. The variation in (stabilizing) axial curvature ∂ xx d around the wave crest (marked by a symbol) significantly increases in magnitude with increasing

Figure 9 .

 9 Figure9. Precipitation of occlusion due to axial viscous diffusion. Cases 1, 2, and 3. Circles/asterisks mark limit points (LP) and absolute instability onset (AI). Solid lines: full model (2.5); dot-dashed lines: no axial viscous diffusion, Jj =Kj=Lj =Mj=0 in (2.5) and (2.11); black curves with symbols: f =fmax; blue lines without symbols: Λ=12.56. The capillary number Ca scales the observed effect. (a,c) Case 1: Ca=0.37. Cross and diamond mark points used in figure10. Panel c represents ratio of interfacial fluid velocity u(dmin) to wave speed c. Grey region marks existence of roll waves u(dmin) > c; (b) case 3: Ca=0.028. In the grey region, the primary flow is absolutely unstable (AI); (d) case 4: Ca=0.012.

Figure 10 .

 10 Figure10. Mechanism delaying occlusion due to axial viscous diffusion. Data correspond to the two solution points marked by a cross and a diamond in panel 9a. Solid lines: diamond in panel 9a, Re=1.47×10 -4 ; dashed lines: cross in in panel 9a, Re=1.0×10 -4 . (a) Surface profiles with zones of negative (grey) and positive (white) Ψxx; (b) normalized differential normal axial viscous stress Ψxx according to (3.2). Axial normal viscous stresses counteract gravity in the wave fronts, tending to steepen the trailing and flatten the leading wave front.

F

  xx is the (dimensionless) axial viscous force acting on the cross section of the liquid film. Its differential dF xx =∂ x F xx dx yields the resulting axial normal viscous force acting on a slice of liquid film. In (3.2), we have normalized dF xx with the differential gravitational force to obtain Ψ xx . Grey and white areas underneath the solid curves in panels 10a and 10b identify regions where the axial normal viscous force is directed counter to (Ψ xx < 0) or in the direction of (Ψ xx > 0) gravity. Based on this, the arrows in panel 10a illustrate the resulting action of axial normal viscous stresses. This action tends to steepen the trailing front of the surface wave (where the residual film and the wave hump are pushed toward one another) and elongate the leading wave front (where residual film and wave hump are pulled away from one another). It counters the gravity-induced compression of the leading wave front and elongation of the trailing wave front and thus weakens the saturation mechanism discussed in section 3.3. The strength of the effect increases with increasing Reynolds number, as evidenced by comparing the two profiles in panel 10b and this explains why the solid and dot-dashed lines in panel 9a diverge with increasing Re. Axial viscous diffusion also greatly affects the wave speed and the streamline pattern in the wave-fixed reference frame. In panel 9c, we have replotted the continuation curves from panel 9a in terms of the fluid velocity at the wave crest u(d min ) normalized by the wave celerity c. When u(d min ) > c, a moving-frame vortex develops in the wave hump.

Figure 11 .

 11 Figure 11. Streamlines in the wave-fixed reference frame for the marked limit points (LP) in panel 9c. Blue lines correspond to liquid film and red lines to core gas flow. (a) Full model (2.5): filled circle in panel 9c. The wave hump contains a vortex; (b) no axial viscous diffusion, i.e. Jj =Kj=Lj =Mj=0 in (2.5): asterisk in panel 9c. No vortex is predicted in the wave hump.

Figure 12 .

 12 Figure 12. Precipitation of occlusion due to inertia. Case 4: low-viscosity silicone oil, Bo=0.993. Solid lines: full model (2.5); dot-dashed lines: no inertia, Si=Fij =Gij=0 in (2.5) and (2.11). (a) Nonlinear travelling-wave solutions. Black lines with symbols: f =fmax; blue lines without symbols: Λ=12.56; (b) dispersion curves of linear spatial growth rate αi: Re=15.4.

Figure 13 .

 13 Figure 13. Precipitation of occlusion due to a counter-current gas flow. Case 4. (a) Travelling-wave continuation at f =fmax. Solid: aerostatic pressure drop, Ψ=1; dashed: counter-current gas flow, Reg=-17.27; solid blue line without symbol: Πµ=0 in (2.5), Reg=-17.27; dot-dashed: Πρ=0 in (2.5), Reg=-17.27; (b) continuation at V l =const and f =fmax. Square: V l /π/R 3 =2.435; diamond: V l /π/R 3 =2.770; (c) wave profiles corresponding to the asterisk and filled circle in panel a; (d) corresponding profiles of gas pressure gradient ∂xpg according to (2.5c) (symbols mark position of wave maximum). Grey/white zones between the dashed profiles in panels c, d demarcate regions of negative/positive ∂xpg. Arrows in panel c indicate action of gas pressure gradient on the liquid film.

Figure 14 .

 14 Figure 14. Streamlines in the wave-fixed reference frame corresponding to the marked limit points (LP) in panel 13a. Within the film, liquid below the vortex flows from right to left. (a) Aerostatic gas pressure drop (filled circle in panel 13a): Ψ=1, Re=18.9, Λ=6.3; (b) counter-current gas flow (asterisk in panel 13a): Reg=-17.27, Re=14.6, Λ=5.7. The gas, flowing from right to left, is constricted by the wave hump, similar to a de Laval nozzle.

Figure 15 .

 15 Figure15. Different routes to occlusion in a spatially evolving film. Computations of experimental runs 20 (panels a,c) and 13 (panels b,d ) in[START_REF] Dao | Experimental study of wave occlusion on falling films in a vertical pipe[END_REF]. (a,b) Travelling wave solutions (right abscissae) at Re=const, and dispersion curves of the linear spatial growth rate αi (left abscissae). Grey zones distinguish linear waves that have no saturated travelling state; (c,d) spatio-temporal computations on an L ⋆ =1 m domain with inlet noise (2.6b); (a,c) experimental run 20: case 2 in table 1, Re=0.0497 (regime of certain occlusion, scenario I); (b,d) experimental run 13: case 3 in table 1, Re=0.258 (regime of conditional occlusion, scenario II). Arrow in panel d marks a coalescence event triggering occlusion.

Figure 16 .

 16 Figure 16. Computations of case 3 (Re0=0.075, Remax=0.39; panel 5b) at lower Re, showing a downstream displacement of the occlusion point w.r.t. panel 15d. Spatio-temporal computations on an L ⋆ =2 m domain with inlet noise (2.6b). (a) Re=0.14: regime of conditional occlusion. Occlusion (through scenario II) occurs much further downstream than in panel 15d; (b) Re=0.07<Re0: regime of impossible occlusion. Surface waves are extremely small.

Figure 17 .

 17 Figure 17. Different occlusion scenarios for case 4 (Re0=1.5, Remax=18.9; panel 6b). Spatio-temporal computations on an L ⋆ =1 m domain with inlet noise (2.6b). (a) Re=19>Remax: regime of certain occlusion (scenario I); (b) Re=5: regime of conditional occlusion (scenario II); (c) Re=1.25<Re0: regime of impossible occlusion. Surface waves are extremely small

Figure 18 .

 18 Figure18. Preventing occlusion in the conditional occlusion regime, Re0<Re<Remax, through coherent inlet forcing. Spatio-temporal simulations from panels 16a and 17b with additional coherent inlet forcing of frequency f and amplitude ǫ1=0.1 (2.6b). All other parameters, including the noise level ǫ2, remain unchanged. (a) Case 3 (Re0=0.075, Remax=0.393; panel 5d): Re=0.14, f =0.5 fmax=0.525; (b) Case 4 (Re0=1.5, Remax=18.9; panel 6b): Re=5, f =0.5 fmax=1.27 Coherent inlet forcing produces a regular train of saturated-amplitude travelling waves, allowing to prevent occlusion versus the natural evolution in panels 16a and 17b.

  Appendix A: Validation of WRIBL model based on transient DNS

Figure 19 Figure 19 .

 1919 Figure19represents comparisons between our WRIBL model (2.5) and DNS using the solver Gerris[START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. All computations were performed for case 4 (low-viscosity silicone oil), which is the most challenging. The solid and dashed curves in panel 19a represent travelling-wave solutions obtained with our model by continuing the liquid volume V l at fixed wavelength Λ=5.4 for two different gas flow situations. The dashed line corresponds to an aerostatic gas pressure drop Ψ=1, and the solid line to a counter-current gas flow at Re g =-17.3. We compare these data with two punctual DNS at V l /π/R 3 =2.35,

Figure 20 .

 20 Figure 20. Travelling waves under the effect of a counter-current gas flow: case 4, f =fmax. (a) New representation of the curves from panel 13a in terms of the normalized pressure gradient Ψ=∆pgFr 2 . Solid: aerostatic pressure drop, Ψ=1; dashed: counter-current gas flow, Reg=-17.27. Solutions between circle and asterisk (red segment) are stable. Filled circles correspond to limit points (LP) in panel 13a; (b) variation of the gas Reynolds number Reg at Re=15.72. The intended continuation path is indicated by an arrow in panel a.

Figure 21 .

 21 Figure 21. Version of the computation from panel 4a with axial viscous diffusion deactivated:Jj =Kj=Lj=Mj =0 in (2.5). The film surface is represented at the same time as in panel 4a.

Figure 22 Figure 22 .

 2222 Figure 22 represents results of a transient periodic computation based on our WRIBL model (2.5) with the additional force term Π θ (4.1). Liquid properties and tube radius correspond to case 4, which, due to the relevance of inertia, is the most challenging. Further, Λ=5.40, V l /π/R 3 =2.80, d plug =0.01R, and λ=1. Panel 22a represents the time evolution of the film surface during formation of a pseudo-plug, and panel 22b represents the fully-developed pseudo-plug. To better guide the eye, we have reproduced and shifted the periodic surface profile in these panels. In contrast to the crude numerical core radius limitation employed in our open-domain computations, our improved representation of liquid plugs ensures volume conservation, and thus the pseudo-plug in figure 22 does not reopen. On the downside, this new representation requires a finer spatial numerical resolution. Thus, open-domain computations, such as the one in panel 17a, become more costly (about four times for case 4). This can

Table

							-3	1.4×10 -6
	2	1260.7	1.049	64.0	4.765	0.30	0.95×10 -3	1.7×10 -5
	3	1223.9	0.167	65.0	3.175	3.53	0.98×10 -3	1.1×10 -4
	4	900	0.0045	20.0	1.5	121.4	1.3×10 -3	4.0×10 -3
	Case Re max (-)	Fr (-)	We (-)	Bo (-)	La (-)	Ca (-)	figures
	1	1.91×10 -4 2.7×10 -3 8.4×10 -5	11.06 6.3×10 -4	0.37	3-5,9,11
	2	0.038	0.040	7.1×10 -3	4.39	0.35	0.142	5,15
	3	0.39	0.063	7.3×10 -3	1.86	9.1	0.028	5,9,15,16,18
	4	18.9	0.42	0.18	0.993	1.3×10 3	0.012	6,9,12-14,17-19

  Figure3. Spatio-temporal computations of the experiments in panels 3a, 3b, and 3c of[START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], which were performed in an L ⋆ =400 mm tube. Parameters correspond to case 1 in table 1. Computations were performed with our model (2.5) using inlet/outlet conditions and the noisy inlet perturbation (2.6), where we have set ǫ1=0 and adjusted ǫ2 to match the experimental wave growth.

	1
	5
	r
	12	16	20	24	28	32	36	40	44	48	52	56	60	64	68	72	76	80
									x

  , but their model did not account for axial viscous diffusion. Our computations in figures 3 and 4 provide further validation, by taking into account this effect, which greatly improves agreement of our model predictions with the experiments of[START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF], as a result of the high viscosity of the working liquid. For example, when running the computation in panel 4a without axial viscous diffusion, we did not observe any occlusion over the entire experimental tube length (see figure 21 in appendix D), which is in qualitative contradiction to the experiments. Also, we will show in subsection 3.4 that neglecting axial viscous diffusion increases Re max by 74% for this liquid (panel 9a), moving it far beyond the experimental occlusion limit, as quantified in terms of Re.

  These waves are generated by a combination of the Plateau-Rayleigh instability, gravitydriven advection, and the Kapitza instability.Using the low-dimensional model introduced in Dietze & Ruyer-Quil (2015), which had been applied only to gravity-free films, we have performed spatio-temporal computations of spatially-evolving falling films and computed travelling-wave solutions through numerical continuation. Our model extends upon the earlier works of[START_REF] Trifonov | Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes[END_REF],[START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF][START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 2. air-driven flow[END_REF],[START_REF] Zhou | Viscoelastic liquid film flowing down a flexible tube[END_REF],[START_REF] Liu | Stability of viscous film flow coating the interior of a vertical tube with a porous wall[END_REF], and Ding et al.

	(2019) by accounting for axial viscous diffusion and inertia. These effects are needed to
	accurately predict occlusion for the real liquids studied here (cases 1-4 in table 1). We
	have validated our model computations with DNS (appendix A and panel 5d).
	Travelling-wave solutions are lost when the liquid Reynolds number Re is increased
	beyond a limit point (LP), low-frequency/long waves being more dangerous (smaller
	Re LP ) than high-frequency/short waves (greater Re LP ). By numerically reproducing the
	occlusion experiments of Camassa et al. (2014), we have proved that surface waves causing
	occlusion systematically lie beyond this limit of travelling-wave solutions, thus validating
	the conjecture formulated by these authors.
	Depending on how occluding waves emerge in a spatially evolving film, either directly
	from linear wave selection (scenario I), or subsequent secondary instability and wave
	coarsening (scenario II), we have distinguished two possible occlusion scenarios. By
	reproducing several occlusion experiments with our own spatio-temporal computations,

we have shown that occlusion scenario I (panels 4a and 15c, movie M2) applies to the experiments of

[START_REF] Camassa | Viscous film-flow coating the interior of a vertical tube. part 1. gravity-driven flow[END_REF]

and run 20 in

Acknowledgements

This work was supported by the ANR wavyFILM project, grant ANR-15-CE06-0016-01 of the French Agence Nationale de la Recherche. Declaration of Interests. The authors report no conflict of interest.

produce an occlusion of the tube by the liquid film, as evidenced by the divergence of the minimal core radius d min . This confirms that occlusion occurs beyond the limit of travelling-wave solutions and that our model is able to capture the associated dynamics.

Appendix B: Numerical procedure for travelling-wave DNS

The DNS results in panel 5d were obtained by solving the full governing equations (2.1)-(2.3) in the limit Π µ =Π ρ =0 with a pseudo-spectral approach. For this, the liquid velocity and pressure fields are projected on modified Chebyshev polynomials:

where X=2(R-r)/(R-d)-1. The radial velocity v l is obtained by integrating (2.1b), i.e.

The functions φ i (X) and ψ i (X) are linear combinations of Chebyshev polynomials of the first kind T i :

so that φ i (X = -1) = 0, and thus the boundary conditions (2.2) are fulfilled.

Writing the Navier-Stokes equations (2.1a) on the Gauss-Lobato points X i =-cos(πi/n) for i 1, and introducing the wave-fixed coordinate ξ = x -c t, yields 2(n -1) relations:

where D ξ =d/dξ, C j = X -1 X(φ j ) ′ dX, and I j = X -1 φ j dX. Further, we have: (4.4) and thus the free-surface dynamic conditions (2.3b)-(2.3c) and the integral continuity equation (2.5a) are recovered in the limit η → 0. We finally obtain the linear equation system:

A dU dξ = B(U; η), (4.5)