Deformation and dewetting of liquid films under gas jets

We study the deformation and dewetting of liquid films under impinging gas jets using experimental, analytical and numerical techniques. We first derive a reduced-order model (a thin-film equation) based on the long-wave assumption and on appropriate decoupling of the gas problem from that for the liquid. To model wettability we include a disjoining pressure. The model not only provides insight into relevant flow regimes, but also is used to guide the more computationally prohibitive direct numerical simulations (DNS) of the full governing equations, performed using two different approaches, the Computational Fluid Dynamics (CFD) package in COMSOL and the volume-of-fluid Gerris package. We find surprisingly that the model produces good agreement with DNS even for flow conditions that are well beyond its theoretical range of validity. We find that for thicker films dewetting is dominated by the gas jet and the dry spot heals over time when the gas flow is switched off. For relatively thin films, the gas jet is important for initiating dewetting, which is then dominated by the receding contact line motion, and the dry spot remains dry after the gas flow is switched off. There is also a range of intermediate thicknesses where the gas-induced dewetted spot may remain dry or heal over time, depending on how strong the gas flow that produced the spot was. The observed behaviours are explained using a bifurcation diagram of steady-state solutions in the absence of the gas flow. We additionally compare the computational results with experiments and find good agreement.


Introduction
The process of impingement of a gas jet onto a liquid layer is important in numerous industrial applications. For example, it is used in steel production in the basic oxygen furnace process (e.g. Turkdogan 1996;Hwang & Irons 2012), in coating applications in the gas-jet wiping process (e.g. Thornton & Graff 1976;Lacanette et al. 2006) and in immersion lithography to remove water from a photoresist coated wafer (e.g. Berendsen et al. 2012Berendsen et al. , 2013. A closely related process is that of impingement of a gas plasma jet (instead of simply a gas jet) onto a layer of a liquid which appears, for example, in the arc welding process (e.g. Berghmans 1972), in medical applications such as wound healing and skin treatment (e.g. Tian & Kushner 2014;Verlackt et al. 2018) and in environmental applications such as water treatment and disinfection (e.g. Foster 2017).
The impact of gas jets onto layers of liquids has been previously studied mainly for the case when the layer of the liquid is relatively thick. A gas jet impinging onto a liquid layer exerts normal and tangential stresses on its surface, which result in its deformation creating a cavity and flow inside the liquid. Most of the previous research was focused on analysing the shape of the cavity and its stability. An early experimental study was performed by Banks & Chandrasekhara (1963), who identified three regimes, namely, a steady cavity, an oscillating cavity and splashing. They focused on the analysis of steady cavities and suggested scaling approaches to establish a relation between the impact of the jet and the depth of the cavity. Turkdogan (1966) carried out the Banks & Chandrasekhara (1963) experiments with liquids of different densities but focused on the effects of the gas nozzle diameter and the nozzle distance from liquid surface. Cheslak et al. (1969) performed an analysis similar to Banks & Chandrasekhara (1963) and concluded that the occurrence of splashing or a smooth cavity depends on the jet velocity, while the viscosity of the liquid and surface tension are less important. Molloy (1970) studied not only the effect of the gas jet on the cavity, but also the effect of the liquid properties. Previous analytical work investigating the shapes of steady cavities has been mainly based on a conformal mapping approach, in which the flow in the liquid is neglected and the system is assumed to be two dimensional, although both of the assumptions are clearly not valid in practice. The first analytical work using a conformal mapping method was done by Olmstead & Raynor (1964), who studied the cavity shape at relatively small gas velocities in the case of small cavity. Vanden-Broeck (1981) used a similar approach but solved the problem using a different numerical procedure, which allowed analysis of the system for larger gas velocities. A more recent analytical work based on a conformal mapping approach was done by He & Belmonte (2010), who analysed the cavity shape without requiring it to be small. Mordasov et al. (2016) employed the balance equations for forces at the gas-liquid interface and not the balance equation for pressure as was used in most previous studies and obtained good agreement with experiments. Despite previous analytical approaches, detailed understanding of the cavity instability mechanisms is still missing. More recent work on gas jets impinging onto liquids has been mainly focused on experimental and DNS investigations (see e.g. Nguyen & Evans 2006;Solórzano-López et al. 2011;Liu et al. 2015;Muñoz-Esparza et al. 2012;Adib et al. 2018) There has been less investigation for the case when the layer of the liquid is relatively thin. In such a case, if the gas jet flow is sufficiently strong, the film ruptures and dewetting is initiated. Gas-jet induced dewetting of thin liquid films was first considered by Berendsen et al. (2012Berendsen et al. ( , 2013 both experimentally and using modelling, via a reducedorder thin-film equation. In Berendsen et al. (2012), the authors focused on analysing the liquid film rupture times and the influence of surfactants and found good agreement between and experimental and modelling results. In Berendsen et al. (2013), the authors additionally analysed the effect of the movement of the gas jet.
In the present study, we present a comprehensive theoretical and experimental study of the deformation and dewetting of (thin and moderately thin) liquid films in a cylindrical beaker under the influence of an impinging gas jet that is generated by maintaining a constant gas flow rate from a stationary cylindrical tube. To obtain initial insight into relevant flow regimes and timescales of the system, we use a systematically derived thinfilm equation that in the axisymmetric case coincides with the model of Berendsen et al. (2012). The equation is obtained under the long-wave assumption and by decoupling the problem for the gas from that for the liquid, under the so-called quasi-static assumption. This involves modelling the gas-liquid interface as a solid wall for the gas problem, which is valid when the typical velocity in the gas is much larger than that in the liquid, see e.g. Tuck (1975) and also more recent work by Tseluiko & Kalliadasis (2011) and Vellingiri et al. (2015). Unlike Berendsen et al. (2012), we use an iterative procedure for the computation of the normal and tangential stresses exerted by the gas onto the gasliquid interface which allows us to obtain more accurate results. The thin-film equation is used to guide the more computationally expensive DNS for the full coupled system of the governing equations for our setup, which we perform using two different packages, the CFD package in COMSOL (see e.g. Pryor 2011) with a moving mesh interface and the volume-of-fluid Gerris package (see e.g. Popinet 2009). These two packages offer distinct advantages and features improving our understanding of the system. DNS are used to estimate the range of validity of the reduced-order model and allow us, on the one hand, to access regimes that would be inaccessible with the reduced-order model and, on the other hand, to analyse flow characteristics that would be very difficult to image in the experiments reliably (e.g. the velocity field). In addition to studying dewetting, we also analyse the post-dewetting dynamics, when the flow of the gas is switched off. An insight into the expected behaviours for various parameter values is provided by a bifurcation diagram of steady state solutions for the system in the absence of the gas flow.
The manuscript is organised as follows: In § 1 we describe our experimental setup. In § 2 we present the governing equations for the system and derive a reduced order model. Next, § 3 explains our computational framework. In § 4 we present and discuss the results. Finally, in § 5 we give our conclusions.

Experimental setup
A schematic representation of the experimental setup is shown in figure 1. We consider a layer of a liquid in a transparent cylindrical beaker and we study the deformation of the surface of the liquid under the influence of an impinging gas jet. The beaker is 6 cm in height and 3 cm in diameter and it is placed in a transparent square tank to minimise the distortion of the image. The liquid used in the experiments is water at room temperature. The gas jet is generated by maintaining a gas flow at a constant rate from a stationary cylindrical tube (nozzle) with its axis coinciding with the axis of the beaker. The inner diameter of the nozzle is 1.6 mm. The nozzle is connected to a compressed gas tank, and the flow rate is controlled by a mass flow controller (MKS, PR4000B). The gas used in most of the experiments is air at room temperature. The nozzle is fixed with a clamp which can be moved to adjust the distance from the nozzle to the surface of the liquid. We typically consider a distance of 5 mm. The position can be read from a movable calibrated scale.
A high-speed camera (Photon Fastcam, M2.1) coupled with a long-distance lens (Infinity, KC) is placed on one side of the square tank in order to record images of the deformed liquid layer. The camera is connected to a computer to enable gathering of the data for analysis. A light source (Kern Dual Fiber Unit LED) is placed on the opposite side of the square tank to provide illumination for the images. The camera is fixed on an adjustable x-y-z stage allowing us to modify the camera's position properly and capture images in the beaker at different places. In particular, the initial position of the camera is adjusted by placing a graticule at the centre of an empty beaker and moving the camera along the stage until the image of the graticule comes into focus. This also allows measuring the size of the interrogation window. An example of an image of the graticule is included in the supplementary material. The recorded images were analysed with the software pack-    Figure 3. Schematic representation of a gas jet impinging onto the surface of a liquid in a cylindrical beaker.
age ImageJ (e.g. Abràmoff et al. 2004). Examples of processed recorded images (with increased contrast) are shown in figure 2 for a relatively thin water film (of undisturbed thickness 0.2 mm) in panel (a) and a relatively thick water film (of undisturbed thickness 1 mm) in panel (b). The corresponding raw images are included in the supplementary material. In panels (a) and (b), the water films were deformed by air jets of flow rates 0.15 and 0.5 slpm (standard litres per minute), respectively. The interrogation window for these image was 2.6 mm × 2.6 mm. A deformation of the gas-liquid interface can be clearly seen in both cases. The gas-liquid interface is constructed by curve fitting in the software package Matlab.

Problem statement and full governing equations
A schematic representation of the model system is shown in figure 3. We denote the radius of the beaker by R b and its height by H. The thickness of the undisturbed liquid layer is denoted by h 0 . The inner and outer radii of the nozzle are denoted by R i and R o , respectively, and the distance between the nozzle and the undisturbed gas-liquid interface is denoted by h 1 . We introduce cylindrical polar coordinates (R, ϕ, z) with the z-axis pointing upwards along the axis of the beaker in the direction opposite to gravity g, and with the bottom of the beaker coinciding with the z = 0 plane. The deformed gas-liquid interface is denoted by Σ and is given by the equation f (R, ϕ, z, t) = 0. In the simplest case, the interface is given by the graph of a function, z = h(R, ϕ, t), and hence f = h − z. We assume that the liquid and the gas are of the same constant temperature. As the gas jet strikes the liquid surface at the centre, the surface of the liquid deforms and a cavity appears. Motion in the liquid, in the form of eddies is also generated. The eddies affect the mass transfer and mixing in the liquid. In order to describe such a system, the Navier-Stokes and continuity equations are used and the corresponding boundary conditions must be satisfied.
As the typical velocity in the liquid will be assumed to be small compared to the speed of sound in the liquid, it is appropriate to model the liquid problem by the incompressible Navier-Stokes equations: where ρ l is the liquid density, which is assumed to be constant, u is the velocity field in the liquid, D/Dt = ∂/∂t + u · ∇ is the usual material derivative, and σ l is the viscous stress tensor for the liquid given by where p l is the pressure in the liquid, µ l is the viscosity of the liquid, which is assumed to be constant, δ is the identity tensor and S l = 1 2 (∇u) T + ∇u is the strain-rate tensor. As we will consider relatively high gas velocities (of up to 50 m/s), compressibility effects may become important in the gas. However, for such gas velocities it is still appropriate to assume that the gas flow is laminar. The Navier-Stokes equation in the gas then takes the form where ρ g is the gas density, v is the gas velocity and σ g is the viscous stress tensor for the gas given by where p g is the gas pressure, µ g is the viscosity of the gas, which is assumed to be constant, and S g = 1 2 (∇v) T + ∇v is the strain-rate tensor in the gas. The continuity equation for the gas is In the gas, we need an additional equation of state, and we impose the ideal gas relation p g = R g ρ g Θ, where R g is the gas constant and Θ is the temperature. We impose no-slip and no-penetration conditions both for the liquid and for the gas at the solid boundaries (u = 0 and v = 0) except at the beaker side wall, where nopenetration still applies, but instead of no-slip, we impose the Navier slip condition to allow for the motion of the contact line (see, e.g., Sibley et al. 2012). So for the liquid when R = R b we have wheren w is a unit normal vector to the side wall pointing into the beaker (i.e.n w = −R, whereR is a unit vector pointing in the R direction),t w is a unit tangent vector to the wall, β l N S is the slip coefficient for the liquid. For the gas when where β g N S is the slip coefficient for the gas. Note that in (2.6) and (2.7), there are, in general, two independent tangent directions to the wall, i.e. each of the equations results in two scalar equations, by takingt w =k andt w =φ, wherek andφ are unit vectors pointing in the z and ϕ directions. However, we will later assume axisymmetry, i.e. no dependence on ϕ, and then it will be sufficient to only considert w =k.
In addition, we impose a fixed contact angle condition at the contact line, so that when the interface is given by z = h(R, ϕ, t), we have where θ c is the angle the liquid makes with the wall at the contact line. In a subset of our numerical simulations presented below, at the bottom of the beaker we also impose the Navier slip condition instead of the no-slip condition. This allows us to study dewetting induced by the gas jet using DNS where topological transitions of the gas-liquid interface are allowed. We do this with the volume-of-fluid package Gerris (e.g. Popinet 2009), as will be explained below. The contact angle at the bottom of the beaker will be denoted by θ eq . We also performed DNS in the CFD finite-element package COMSOL (e.g. Pryor 2011), and our implementation allows for mesh movements so that the mesh deformations follow the gas-liquid interface motion. Such an implementation allows us to analyse the deformations of the interface very accurately but forbids topological transitions.
At the gas inlet, when z = h 0 + h 1 and 0 R R i , we impose the fully developed laminar Poiseuille velocity profile: where v max = 2q g /πR 2 i , with q g denoting the imposed gas flow rate. At the gas outlet, when z = H and R o < R < R b , we impose normal flow, and prescribe normal stress, i.e. we requirek · σ g ·k = −p a , where p a is the atmospheric pressure.
Finally, we discuss conditions that must be satisfied at the gas-liquid interface Σ. First, we have the kinematic condition Df Dt = 0, (2.10) where we remind that f is a function such that Σ is given by the equation f (R, ϕ, z, t) = 0. Continuity of velocity must also be satisfied at the interface, u = v, and we must have dynamic balance of stress at Σ:n · σ l −n · σ g = γκn. (2.11) Here,n is the unit normal vector to the interface pointing into the liquid. The term on the right-hand side is due to the Laplace pressure, where γ is the gas-liquid surface tension coefficient (which is assumed to be constant) and κ = ∇ ·n is twice the mean curvature of the interface Σ. Note that to study dewetting induced by the gas jet in a numerical formulation where topological transitions are not allowed, we also include a Derjaguin (or disjoining) pressure in the stress balance condition. This approach is applicable when Σ is a graph of a function, z = h(R, ϕ, t). The stress balance condition then becomeŝ n · σ l −n · σ g = γκn + Π(h)n. (2.12) The disjoining pressure represents an effective interaction between the gas-liquid interface and the liquid-substrate interface. It can be written as Π(h) = −dV ( where the first term results from the long-range attractive forces (with A representing the Hamaker constant) and the second term results from the short-range repulsive forces. The second term prevents the liquid film from breaking down, and instead of this occurring we obtain a very thin precursor film. In practice, where the film thickness is equal to the thickness of the precursor film, we may assume that a dry spot has appeared. At equilibrium, the thickness of the precursor film corresponds to the minimum of the binding potential V (h) and is equal to h eq = (B/A) 1/3 , with the contact angle at the apparent contact line is given by (e.g. Rauscher & Dietrich 2008;Hughes et al. 2015) θ eq = cos −1 1 + V (h eq ) γ . (2.14)

Thin-film model
Solving the full system of governing equations is a computationally expensive task. We therefore aim to simplify the problem by deriving an accurate reduced-order model. Such a model not only provides insight into the fundamental features of the system in an efficient way but also serves as a mechanism to guide the more computationally prohibitive DNS tools towards suitable regimes with a much more informed view of an otherwise vast parameter space. The first step is to decouple the problem for the gas from that for the liquid, which is possible when the typical velocity in the liquid is much smaller than that in the gas. Then, for the gas problem it is appropriate to neglect the motion of the liquid and to use the quasi-static assumption, i.e. it is appropriate to model the interface as a rigid wall and solve the gas problem independently, see e.g. Tuck (1975) who states that such an assumption is appropriate when the typical liquid velocity is less than approximately 4% of the typical gas velocity, which is always the case in our study (see also Tseluiko & Kalliadasis 2011;Vellingiri et al. 2015, for more recent studies, where the quasi-static assumption was used in the modelling of a liquid film sheared by a turbulent gas). The solution of the gas problem can then be used to obtain the stress exerted by the gas onto the gas-liquid interface, which can then be fed into the normal and tangential stress balance conditions at the interface. We denote such a stress by s g so that For the analysis in this section, we first consider the general non-axisymmetric case, and for convenience we use Cartesian coordinates (x, y, z) (so that the z direction remains as before. We non-dimensionalise the equations using h 0 as the length scale, U 0 as the velocity scale (to be specified later), h 0 /U 0 as the time scale and µ l U 0 /h 0 as the scale for pressure and the gas stress, so from now on all the variables will be assumed to be dimensionless. We thus introduce dimensionless variables via the following mappings: where we denote by u, v and w the x, y and z components of the velocity, respectively.

The incompressible Navier-Stokes and continuity equations in the liquid become
where Re and G are the Reynolds and the gravity numbers, respectively, given by The no-slip and no-penetration conditions at the bottom of the beaker become The kinematic condition at the interface, z = h(x, y, t), takes the form The normal stress balance condition takes the form: where Ca is the Capillary number given by N s is the dimensionless normal stress exerted by the gas on the gas-liquid interface which under the quasi-static assumption can be assumed to be a functional of the interface shape, N s = N s [h], and finally, Π(h) is the dimensionless disjoining pressure, given by x andt =t 2 = (0, 1, h y )/ 1 + h 2 y in the tangential stress balance condition, we obtain where T si , i = 1, 2, are thet 1 andt 2 components of the tangential stress exerted by the gas on the gas-liquid interface, T si = s g ·t i , which under the quasi-static assumption can be assumed to be functionals of the interface shape, T si = T si [h]. The tangential stress exerted by the gas on the gas-liquid interface is then expressed as Next, we utilise the so-called thin-film or long-wave approximation, namely, we assume that the undisturbed film thickness, h 0 , is much smaller than the characteristic horizontal length scale over which variations in the film thickness occur, and we introduce the socalled thin-film parameter = h 0 / 1. We now use the following additional rescalings of variables that are standard for the thin-film approximation: To derive the thin-film equation, we consider the asymptotic limit → 0. Then, to keep capillary effects at leading order, we assume that Ca is asymptotically bounded above and below by 3 . To neglect inertia at leading order, we assume that Re 1/ . We also assume that G = O(1/ ), so that gravitational effects may enter at leading order. For the disjoining pressure, it is appropriate to assume that Π = O(1/ ). In addition, for the dimensionless gas stress, we need to assume that N s = O(1/ ) and T si = O(1), i = 1, 2. We then introduce the following rescaled parameters: so that Ca is asymptotically bounded above and below by non-zero constants and G = O(1), and the following rescaled gas normal stress: , and the rescaled disjoining pressure: The problem at leading order becomes: with the no-slip and no-penetration conditions at z = 0 and the tangential and normal stress balances We also have the kinematic condition, which can be rewritten as is the flux vector parallel to the plane z = 0. Then we find that the pressure at leading order is given by 41) and the velocity components at leading order are given by Substituting the leading-order expressions for u and v into the expression for the flux vector (2.40), we find that at leading order Substituting this expression into the kinematic condition (2.39) gives the following evolution equation for the film thickness, the so-called thin-film equation: Scaling back to the dimensionless variables x, y and t, we obtain where the dimensionless leading-order pressure is given by and ∇ = (∂/∂x, ∂/∂y). For convenience, it is possible to eliminate one of the dimensionless parameters by, for example, multiplying the thin-film equation by Ca and appropriately rescaling time. This is equivalent to choosing U 0 = γ/µ l , for which the time scale becomes µ l h 0 /γ and the scale for pressure and the gas stress becomes γ/h 0 . Then the pressure in the thin-film equation takes the form where Bo is the Bond number given by and the dimensionless coefficients in the disjoining pressure are A = A/γh 2 0 and B = B/γh 5 0 . For the validity of the thin-film equation in terms of dimensionless parameters that are independent of U 0 , we must have Bo = O( 2 ) and La 1/ 4 , where La is the Laplace number given by (2.50) The latter condition is needed for inertia to be negligible. We must additionally have The validity of the thin-film equation for the experimental parameter values that we have used is discussed § 4. Finally, going back to cylindrical polar coordinates (R, ϕ, z) and assuming axisymmetry, we obtain the following equation: where T s denotes the gas tangential stress in the R-direction and pressure is given by Note that here R is assumed to be non-dimensionalised using h 0 as the length scale.
To solve the thin-film equation (2.51) numerically, we also need to impose appropriate boundary conditions. Conditions at R = 0 follow from the symmetry assumption: At the side wall, we will assume for simplicity that the contact angle is 90 • , so that For analysing flow patterns in the liquid film, it is also useful to give the R and z velocity components in cylindrical polar coordinates: (2.57)

Computational framework
Complementing the experimental and analytical investigations we also considered two distinct numerical platforms to simulate the target physical system. The two packages (described in more detail in the paragraphs to follow) offer distinct advantages and features that aid our understanding of the flow characteristics. They act not only to bridge the gap between the previous approaches, but also to access regimes that would be inaccessible with a reduced-order model approach on the one hand, as well as easily allow the inspection of quantities in the flow that would be very difficult to image reliably on the other.
First, we implemented the setup in the commercial software platform COMSOL Multiphysics 5.3a. We used the CFD module which is a standard tool to simulating systems that involve complex fluid flow models. A two-dimensional axisymmetric geometry was built using the parameters from the experiments. COMSOL uses an unstructured mesh finite-element approach, which is highly suitable for tracking details near specific target regions of the domain. However, for the present problem, we found it challenging and computationally highly expensive to accurately describe the evolution of the gas-liquid interface and topological transitions, occurring e.g. in the dewetting process, using the built-in level-set and phase-field methods. We thus utilised a more computationally efficient moving-mesh approach in which the gas-liquid interface is modelled as a sharp surface separating the two phases and the mesh deformations follow the deformations of the interface. However, such an implementation is not directly suitable for describing topological transitions such as in the dewetting process. Thus, as discussed in § 2.1, we included the disjoining pressure into the normal stress balance condition to study dewetting which prevents the liquid film from breaking down so that a dry spot is modelled with a very thin precursor film. We should note that this approach still has limitations in modelling dewetting. For example, it is not suitable for describing dewetting on hydrophobic surfaces for which the contact angle is greater than 90 • .
To overcome the limitations of our COMSOL implementation with respect to topological transitions, we also implemented the setup in the open-source package Gerris (e.g. Popinet 2009). Well-known in the interfacial flow community for more than a decade, its strengths lie in the adaptive mesh refinement and parallelisation capabilities that make it an ideal testbed for multi-scale flow problems. The transparent structure of the code allows for careful validation of any in-house implemented extensions, as have been employed here. For example, one particular region of interest in the flow is the near-wall region where dewetting can be considered without the need to introduce a precursor film. The interface-capturing techniques, coupled with well-established contact line models (e.g. Afkhami et al. 2018) and control over any imposed Navier-slip-type conditions, provide an added perspective to the overall investigation. The chosen refinement strategy concentrates on adequately addressing the sensitive regions near the gas nozzle and the walls in contact with the liquid, while adaptive refinement is used to steer degrees of freedom towards any changes in interfacial position, as well as changes in components of the velocity field and vorticity in order to accurately capture non-trivial flow regions in both the liquid and the gas. The runs in the sections to follow have been executed in parallel on local computing facilities, typically amounting to O(10 3 ) CPU hours in Gerris, depending on flow parameters. While the chosen adaptive strategy restricts the number of gridnodes to O(10 5 ), a gain of two orders of magnitude over a uniform grid with the same minimum cell size, the delicate interplay between the gas-liquid coupling requires special measures from a linear algebra and stability viewpoint, leading to relatively small time steps in exchange for mesh-independent results. By contrast, we found that the computations in COMSOL were significantly faster for a similar number of degrees of freedom (we typically used O(10 5 ) triangular elements), with a typical computation for the fully coupled problem taking approximately O(10 2 ) CPU hours. However, as already noted above our implementation in COMSOL has its restrictions for systematically studying dewetting. We have validated both implementations extensively by systematically decreasing the cell size or increasing the number of mesh elements until the convergence in the numerical results was achieved. The results of the computations and comparisons with the experimental and reduced-order model results are presented in the next section.

Results
Throughout this section, we consider the following geometrical parameters in both the mathematical model and in the numerical simulations: the diameter of the beaker is 3 cm and its height is 6 cm; the inner diameter of the gas nozzle is 1.6 mm and its distance from the undisturbed liquid surface is 5 mm, as in the experiments. Also, the liquid is water and the the gas is air at room temperature. The density and viscosity of water at room temperature are ρ l = 1000 kg m −3 and µ l = 8.9 × 10 −4 Pa s, respectively. For the density and viscosity of air, we use ρ g = 1.22 kg m −3 and µ g = 1.81 × 10 −5 Pa s. The surface tension coefficient for the air-water interface is set to γ = 72 × 10 −3 N m −1 . Regarding the results obtained with the thin-film equation (2.51), for the asymptotic validity of the equation Bo 1 is required and also La 1/Bo 2 . It can be verified that for a water film La becomes smaller than 1/Bo 2 if h 0 < 0.226 mm (and then Bo < 0.007). Thus, strictly speaking for the validity of the thin-film model the thickness of the film must be less than ≈ 0.23 mm. As regards gas flow rates suitable for the validity of the thin-film equation, we note that for a film of thickness 0.226 mm and an air jet flowing at the rate 0.2 slpm, the maximum values of the normal and tangential stresses non-dimensionalised with γ/h 0 turn out to be 0.007 and 0.0008, respectively, which is appropriate for the validity of the model. For an air jet flowing at the rate 0.4 slpm, the maximum values of the dimensionless normal and tangential stresses turn out to be 0.028 and 0.0025, respectively, which is also close to the region of the validity of the model. However, it will be shown below that the thin-film equation turns out to produce good agreement with DNS and experiments also for film thicknesses significantly larger than 0.226 mm (of O(1) mm) and for gas flow rates significantly higher than 0.4 slpm (of O(1) slpm). For even thicker water films and higher gas flow rates, inertial effects become important and the derived thin-film equation is not valid. Such parameter regimes can be accessed with the developed DNS framework.

Decoupled gas problem
First, we discuss the decoupled gas problem and explain how the gas normal and tangential stresses exerted on the interface are computed using an iterative procedure. Under the quasi-static assumption, we model the liquid surface as a solid wall obtaining a problem for the gas only, and initially we assume that the interface is flat. The gas flow rapidly develops into a steady state. The resulting gas flow pattern at steady state obtained using COMSOL is shown in figure 4(a) for the gas flow rate q g = 1 slpm, which corresponds to the maximum gas speed of approximately 16 m s −1 . In this section, we assume that the undisturbed interface is located at z = 0. We note that Gerris simulations agree with the COMSOL results. The colour scheme indicates the gas speed in metres per second and the thin white lines show streamlines. We can observe that the gas jet impinges onto the lower wall, and then the gas flows in the direction parallel to the wall radially outwards with its speed decreasing as the radial distance increases. We can also observe that a relatively large recirculation zone (eddy) is generated in the gas, and there is also a small, relatively slow eddy in the bottom-right corner. This solution of the gas problem for the interface modelled as a flat solid wall is used to compute the normal and tangential stresses exerted by the gas jet on the interface at different radial locations. Next, we use these stresses in the thin-film equation (2.51) to solve the problem for the liquid film and therefore obtain the deformation of the interface resulting from these stresses. The thin-film equation is solved numerically in Matlab using finite-difference approximations for the spatial derivatives and Matlab's ode15s solver for stepping in time. For this gas flow rate, the solution evolves into a steady state within a few seconds. An example of a deformed interface computed in this way is shown in figure 4(b). A liquid film of thickness 5 mm was used for illustrative purposes (although we note that the thin-film equation is not expected to be valid for such a thickness). Next, we use this deformed interface as the lower boundary for the gas domain and again assume that it is a solid wall and recompute the solution of the decoupled gas problem. It can be seen in figure 4(b) that the computed solution is qualitatively similar to the one in figure 4(a). We extract the gas stresses from this solutions and use these updated stresses to solve the thin-film equation again to obtain an updated steady-stated interface shape. This procedure is repeated until a converged steady-state interface shape is achieved. Typically, we find that 2-3 iterations are sufficient.
Next, we will analyse in detail how the normal and tangential stresses exerted by the gas onto the interface behave when the gas flow rate varies for the case when the interface is modelled as a flat solid wall. The results are presented in figure 5. Panels (a) and (d) show the normal and tangential stresses for the gas flow rate q g changing from 0.2 (red lines) to 1.8 slpm (blue lines) with the increment of 0.2 slpm. As expected, the stresses grow as the gas flow rate increases. The normal stresses have their maximum values in the centre (at R = 0) and then rapidly decay as the radial distance increases. The tangential stresses vanish at the centre and have their maximum values at a distance slightly away from the centre, at approximately R = 0.1 cm. Then they slowly decay as R increases up to approximately R = 1.2 cm. After this distance, the tangential stresses become negligibly small, and this may be associated with the presence of a slow recirculation zone in the corner, as seen in figure 4. Panels (b) and (e) show the maxima of the normal and tangential stresses (blue lines), respectively, versus the gas flow rate, q g , on the loglog scale and suggest a power law behaviour, i.e. max N s scales approximately as q 2 g and max T s scales approximately as q 1.45 g . Indeed, the red dashed lines in panels (b) and (e) have slopes 2 and 1.45, respectively. It can be observed that the scaling for the normal stresses works well for all the values of the flow rate, whereas for the tangential stresses the scaling works better for larger values of the flow rate. We plotted the rescaled normal and tangential stresses N s /q 2 g and T s /q 1.45 g in panels (c) and (f), respectively. For the normal stresses, we can observe that the curves seem to collapse onto the same universal curve (except for the smallest gas flow rate for which there is a slight deviation, see the red curve). For the tangential stresses, the scaling works very well only for a relatively small range of R values. After the maxima at approximately R = 0.1 cm, the curves start to deviate from each other, but seem to approach a universal curve for larger values of q g . The scaling for the normal stress follows from the fact that the stagnation point pressure (i.e. the pressure or normal stress where the gas jet impinges on the wall at R = 0) is associated with the dynamic pressure at the centre line of the gas jet, which follows from Bernoulli's theorem (assuming incompressible flow), see e.g. Cheslak et al. (1969); Clancy (2006). The scaling then follows from the fact that the dynamic pressure is proportional to the square of the jet velocity, which in turn is proportional to the flow rate q g . The reason for the apparent scaling for the tangential stresses is not immediately obvious from the governing equations and is left as a topic for future investigation. We can conclude that these scalings may be utilised for larger values of the gas flow rate (particularly for thin films where the interface is not deformed much), but for smaller gas flow rates these scalings do not work well and the stresses need to be recomputed for each value of q g (this particularly applies to tangential stresses).

Comparison of the thin-film model with DNS and experiments
In this section, we present results for a film of thickness 0.5 mm. We start with the gas flow rate of 0.2 slpm at which the film deforms but does not rupture. The results are given in figure 6. Panels (a) and (b) show the normal and tangential stresses, respectively, exerted on the interface. These are computed using the iterative procedure described above. The blue solid lines correspond to the flat interface, the black dotted lines correspond to the curved interface, obtained by solving the thin-film equation (2.51), after the first iteration, and the red dashed lines are obtained for the curved interface after the second iteration. Step-1 curved surface Step-2 curved surface It can be observed that two iterations are sufficient in this case and convergence is achieved. There is only a minor difference in the results for the normal stresses. However, there is a noticeable difference between the tangential stresses computed for flat and curved interfaces. The resulting interface deformations for this gas flow rate are shown in panel (c) after the computational time t = 5 s, although to converge to a steady state approximately 2-3 seconds was found to be sufficient. In panel (c), we compare the thinfilm results computed with the gas stresses obtained using the different stages of the iterative procedure (the brown thick solid and black thin solid lines) with the Gerris and COMSOL results for the fully coupled gas-liquid model (the green dash-dotted and blued dashed lines, respectively). An experimental result is also shown (the red dotted line).
There is a slight difference in the region near R = 0 between the Gerris and COMSOL results for the fully coupled model and the thin-film result when the gas stress were computed by assuming that the interface was flat. This difference is nearly eliminated after the iterative procedure, and we conclude that the thin-film equation performs well even for a film thickness well beyond the expected range of validity of the equation. Panel (d) shows the streamlines in the liquid film obtained using the thin-film model (blue solid lines) and COMSOL (red dotted lines). We can observe that there is a relatively large eddy close to the axis of symmetry at R = 0 and there is a smaller and slower eddy near Step 1 thin-film result Step 2 thin-film result Step 3 thin-film result COMSOL result Figure 7. Numerical results for a gas jet of flow rate qg = 0.37 slpm impinging onto a liquid film of thickness h0 = 0.5 mm. Shown are the resulting interface deformations at steady state obtained using the thin-film equation (2.51) with the iterative procedure and COMSOL for the fully coupled gas-liquid model, as indicated in the legend. The numerical formulations include the disjoining pressure (2.13) giving the equilibrium contact angle θeq = 15 • and the precursor thickness heq = 0.03 mm.
the side wall of the beaker. This second eddy appears due to small and slow eddy in the gas in the corner between the gas-liquid interface and the side wall, see figure 4. We note that the Gerris results for the streamlines are in qualitative agreement with this.
The importance of the iterative procedure in computing gas stresses becomes more apparent when the gas flow rate is close to the value that leads to dewetting of the liquid from the bottom of the beaker. An example is shown in figure 7 for gas flow rate q g = 0.37 slpm. To account for possible dewetting, the thin-film and COMSOL numerical simulations here include the disjoining pressure (2.13) with A = 2.9609 × 10 −11 and B = 8.2659 × 10 −25 giving the contact angle θ eq = 15 • and the precursor thickness h eq = 0.03 mm. It can be seen that in the COMSOL simulation for the fully coupled gas-liquid model a dry spot for radius of approximately 0.36 cm appears in the centre (blue dashed line). However, for the first step of the iterative procedure for the thinfilm model we find that no dry spot appears. This is then suitably accounted for by the next step of the iterative procedure. Indeed, when we recompute the stresses for the resulting deformed interface and use them in the thin-film equation again, we observe that dewetting is initiated and a dry spot appears of radius that is in good agreement with the COMSOL result (see the black thin line). The next iteration on the gas stresses and the thin-film equation leads to the result that is indistinguishable from the one in the previous iteration.

Dewetting
Before discussing dewetting induced by gas jets, it is useful to analyse the linear stability of flat-film equilibrium solutions in the absence of the gas jet as well as to construct bifurcation diagrams of various non-equilibrium solutions. It will be shown later that many of the observed behaviours in gas-jet-induced dewetting can be explained by such an analysis.

Linear stability analysis and equilibrium solutions
We will utilise the thin-film equation (2.51) for our analysis. We consider a liquid film in the absence of a gas jet and analyse steady-state solutions of the thin-film equation, which takes the form Due to the destabilising effect of the long-range attractive forces, a sufficiently thin film must become linearly unstable. In (4.1) distances have been non-dimensionalised with the undisturbed film thickness, h 0 , making the undisturbed dimensionless film thickness is 1, so varying the dimensional film thickness is equivalent to varying the dimensionless parameters (and the dimensionless domain size). Thus, we consider the dimensionless equation (4.1) and perform the linear stability analysis of the solution h ≡ 1 to find out the linear instability conditions in terms of the dimensionless parameters. Note that the linear stability analysis for a similar equation but with the disjoining pressure containing only a destabilising term was performed by Witelski & Bernoff (2000). Substituting h(R, t) = 1+h 1 (R, t) into (4.1) and assuming that |h 1 (R, t)| 1, we obtain the following linearised equation: The associated boundary conditions are It can be concluded from Witelski & Bernoff (2000) that the eigenvalues of the operator L are given by (4.5) where Λ n 's are the eigenvalues of the Helmholtz problem with h R = 0 at R = 0 and at R = R b . The eigenfunctions of L are the corresponding eigenfunctions of the Helmholtz problem. The eigenvalues of the Helmholtz problem are real, and it is enough to consider the positive ones. Assuming that the smallest positive eigenvalue is Λ 1 , the condition for linear instability becomes λ 1 > 0, or, equivalently, The eigenvalues of the Helmholtz problem are solutions of the equation J 1 (ΛR b ) = 0, where J 1 is the first-order Bessel function of the first kind. Denoting the smallest positive root of J 1 by x 1 ≈ 3.832, we then find that Λ 1 = x 1 R b . The corresponding eigenfunction is J 0 (Λ 1 R), where J 0 is the zeroth-order Bessel function of the first kind. The linear instability condition becomes (4.8) In terms of the dimensional parameters, this condition can be rewritten as For a given equilibrium precursor thickness, it can be verified that the inequality (4.9) has solutions if the contact angle is sufficiently large. For example, taking h eq = 0.03 mm, we find that the liquid film may become linearly unstable if θ 1.21 • , and taking h eq = 0.015 mm, we find that the liquid film may become linearly unstable if θ 0.61 • .
If the latter condition is satisfied, we find that there exist a range of film thicknesses, h 0 ∈ (h a , h b ), for which the flat film solution is linearly unstable. For example, for h eq = 0.03 mm and θ eq = 30 • , we find that h a ≈ 0.03822516 mm and h b ≈ 0.27959481 mm, and for h eq = 0.015 mm and θ eq = 30 • , we find that h a ≈ 0.01911091 mm and h b ≈ 0.19778522 mm. (Overall, h a and h b become smaller as h eq decreases for a given value of θ eq .) According to the bifurcation theory, for film thicknesses close to the values h a and h b at which linear stability of the uniform-thickness solution changes, there exist non-uniform solutions that bifurcate from the uniform solutions. (There of course exist other branches of solutions that bifurcate from the points where more modes become unstable. Such solutions are, however, less relevant to the present study and we do not discuss them here.) The nature of the bifurcations can be analysed by utilising a weakly nonlinear analysis to obtain the evolution equation for the amplitude of the unstable mode J 0 (Λ 1 R), see e.g. Witelski & Bernoff (2000) for a similar analysis. It turns out that the bifurcations are transcritical, and this is demonstrated e.g. in figure 8(a) for h eq = 0.03 mm and θ eq = 30 • , which shows the bifurcation diagram of the non-uniform solutions when h 0 is used as the bifurcation parameter. We use the dimensionless quantity [h(0) − h(R b )]/h 0 as a measure of the interfacial shape distortion of the solutions. The uniform solutions then correspond to the horizontal line with the measure equal to zero (see the black line). The solid lines show stable solutions and the dashed lines show unstable solutions. The black filled circles show the primary bifurcation points. We find that these points are connected with each other by branches of non-uniform solutions, and, in fact, we find that from each of the primary bifurcation points there emerge two types of solutions. Namely, we obtain solution which have a maximum at R = 0 (the corresponding branch is above the horizontal line at zero and is shown in red). For h 0 = h a , we find that this branch initially goes to the left and is initially unstable up to the first turning point (the turning points are shown with empty circles). For h 0 = h b , we find that this branch also initially goes to the left and is initially stable. We also obtain solution which have a minimum at R = 0 (the corresponding branch is below the horizontal line at zero and is shown with blue colour). Solutions of this type are more relevant to the study of the deformation of liquid films under gas jets, and we will, therefore, discuss them in more detail. For h 0 = h a , we find that this branch initially goes to the right and is initially stable for a very small range of values of h 0 (up to the turning point at h 0 = h t1 ≈ 0.03822585 mm). After the turning point at h 0 = h t1 the branch becomes unstable up to the second turning point at h 0 = h t2 ≈ 0.03355743 mm. Then, the solutions become stable up to the next turning point at h 0 = h t3 ≈ 0.9016 mm. The next part of the bifurcation curve is unstable and goes to the left terminating at the primary bifurcation point, h 0 = h b . The insets in the figure zoom into the regions around the primary bifurcation points and confirm the transcritical nature of the bifurcations.
We conclude that for h 0 > h t3 stable non-uniform solutions having a minimum in the centre do not exist. Thus, if we consider a film of an average thickness h 0 > h t3 that was initially deformed by a gas jet with the gas flow switched off after some time, we expect the dry spot to heal with time with the liquid film eventually levelling up and returning to a uniform-thickness state. Note that the healing process of liquid films has been analysed in detail by Dijksman et al. (2015), Bostwick et al. (2017), Zheng et al. (2018), and in Zheng et al. (2018) in particular it has been shown that there exists selfsimilarity in the healing process and the solutions that govern the healing process are self-similar solutions of the second kind.
For h 0 ∈ (h b , h t3 ), in addition to stable uniform-thickness solutions, there exist stable nonuniform-thickness solutions which have a form of dewetted liquid films with a dry spot in the centre. (As mentioned above, we ignore solutions with a maximum at R = 0 as these are not relevant to our study.) These solutions are separated from the uniform-thickness solutions by an unstable part of the branch of nonuniform solutions. The solutions of this unstable part of the branch also have a form of dewetted liquid films with a dry spot in the centre for h 0 0.38 mm but with the smaller radii of the dry spots compared to the dry spots for the corresponding stable solutions. For h 0 0.38 mm the solutions of the unstable part have a localised minimum in the centre, with the thickness greater than the precursor film thickness, so that a dry spot does not appear. The stable and unstable solutions for h 0 = 0.8 mm are shown in figure 9 with the thick (blue) solid and dashed lines, respectively. As regards a liquid film deformed by a gas jet for such film thicknesses, we expect that if the gas flow rate is not strong enough, the liquid film may initially dewet in the centre but then return to the uniform state after the gas jet is switched off. However, if the gas jet is strong enough to push the liquid film profile beyond the unstable non-uniform solution, we expect that the liquid film remains dewetted in the centre even after the gas flow is switched off. For h 0 ∈ (h t1 , h b ), the uniform thickness solution is linearly unstable and we therefore expect that for any gas flow rate the liquid film will dewet and remain dewetted even after the gas flow is switched off. We are not interested in the behaviour of very thin liquid films of thicknesses close to or smaller than the precursor-film thickness and thus we do not discuss the expected behaviour for h 0 < h t1 .
In figure 8(b), we show the bifurcation diagram of the non-uniform solutions for the same equilibrium contact angle (θ eq = 30 • ) as in figure 8(a) but for a twice as thin equilibrium precursor thickness, h eq = 0.015 mm. We notice that the bifurcation diagram agrees well with the one for h eq = 0.03 mm, particularly for larger values of h 0 . One qualitative difference is the appearance of two additional turning points on the branch of non-uniform solutions with the minimum in the centre, at h 0 = h t4 ≈ 0.262 and at h 0 = h t5 ≈ 0.263. This can be clearly seen in the inset zooming into the region around these turning points. This implies that there is a small range of the film thicknesses, h 0 ∈ (h t4 , h t5 ), where there exist additional stable non-uniform-thickness solutions that have a localised minimum in the centre (without a dry spot). To confirm qualitative and quantitative similarity of the results for h eq = 0.015 mm with the results for h eq = 0.03 mm for larger values of h 0 , we show in figure 9 the stable and unstable solutions for h 0 = 0.8 mm when h eq = 0.015 mm with the thin (red) solid and dashed lines, respectively, in addition to the corresponding solutions for h eq = 0.03 mm. We indeed observe that the results for h eq = 0.03 mm and h eq = 0.015 mm agree very well.

Dewetting induced by gas jets
We first consider a liquid film of thickness h 0 = 0.2 mm. For such a thin film, we used h eq = 0.01 mm for the COMSOL and thin-film computations. To compare with experiments, we measured the equilibrium contact angle for the bottom of the beaker made of an acrylic polymer. This was done by placing a drop of water onto a plate made of the acrylic polymer and measuring the resulting angle using the Drop Shape Analyser DSA100 by KRÜSS. The experimental results showing the evolution of the contact angle θ eq over time are given in figure 10 and indicate that θ eq ≈ 30 • (±2 • ). Using the analysis from the previous section, it can be shown that for such values of θ eq and h eq , a water film of thickness h 0 = 0.2 mm is linearly stable, but this value is close the value h b ≈ 0.1615 mm below which the liquid film becomes linearly unstable. This indicates that a gas jet of a relatively weak flow rate can destabilise the liquid film so that it would dewet leaving a dry spot in the centre of the beaker. This is confirmed by our numerical simulations using all the three different approaches (Gerris, COMSOL and the thin-film model). Figure 11(a) shows dewetted liquid film profiles after t = 1 s for the gas flow rates q g = 0.2, 0.4 and 0.8 slpm using DNS in COMSOL. It can be observed that the profiles converge to the same steady state independent of the gas flow rate. We note that the simulations with the thin-film model produce the profiles which are in excellent agreement with the COMSOL results. The Gerris simulations also agree very well with these results (with a very small difference in the amplitude for the final film profile due to the fact that for the Gerris simulations there is no precursor film and thus the volume of the liquid in the final wetted region is slightly bigger). We, therefore, do not show the thin-film and Gerris results in figure 11(a) (detailed comparisons are discussed below for thicker water films). An experimental result for a dewetted water film of thickness 0.2 mm induced by a gas jet of flow rate 0.4 slpm is given in figure 11(b) which shows a top view for the final profile. The radius of the dry spot turns out to be R d ≈ 1.13 mm, which agrees well with the numerical simulations in which the radius is approximately R d ≈ 1.12 mm in the Gerris calculation and 1.17 mm in COMSOL. In this case, we can conclude that the gas jet is important for initialising dewetting, but dewetting itself is dominated by the receding contact line motion until the equilibrium contact angle is reached. The gas jet does not affect the final steady-state profile. This can be explained by the fact that for the equilibrium solution the wetted region is such that the gas normal and tangential stresses are negligibly small, see figures 5(a,d). Numerical results obtained in COMSOL for the final steady-state profiles for the case of a smaller contact angle and the same film thickness are shown in figure 12 for the gas flow rates q g = 0.2, 0.8, 1.4 and 2 slpm. In this case, the radius of the dry spot for the steady-state solutions is approximately 1 mm, and the influence of the gas jet, although weak, becomes noticeable. As expected, the stronger the gas flow the larger the radius of the dry spot becomes. By looking at the gas normal and tangential stresses in figures 5(a,d), we can notice again that normal stresses are negligible in the wetted area, but now there are small but non-negligible tangential stresses when R ≈ 1 mm, which push the liquid away from the centre. Now we consider a thicker water film with h 0 = 0.5 mm and we assume that θ eq = 30 • and h eq = 0.03 mm. Figure 13(a) shows liquid film profiles after t = 1 s for the gas flow rates q g = 0.2, 0.3, 0.4 and 0.5 slpm (see the black, green, red and blue lines, respectively). Note that at t = 1 s steady state profiles are reached. This can be confirmed in figures 13(b,c) showing the evolutions of the minimum and maximum values of the show reasonably good agreement with the Gerris and COMSOL results, but do not feature oscillations. This may be due to the fact that for such a thickness of the film inertial effects become important and cannot be neglected to accurately describe the evolution of the liquid film. We can observe that dewetting is initiated for a gas flow rate between 0.3 and 0.4 slpm, and, as before, the final dewetted profiles do not depend much on the strength of the gas jet. This also agrees with our experimental observations. Next, we consider even thicker water films with h 0 = 0.8 mm and h 0 = 1 mm, see figures 14 and 15, respectively. As before, we assume that θ eq = 30 • and h eq = 0.03 mm. Panels gas flow switched on and the dashed parts correspond to the gas flow switched off. We can observe that dewetting is initiated for a gas flow rate between 0.5 and 0.525 slpm. However, for q g = 0.525 slpm we find out that after the gas flow is switched off, the dry spot in the centre heals and the liquid films returns to the uniform thickness state. This agrees with the theoretical analysis of § 4.3.1, where we predicted (using the thinfilm equation) the coexistence of stable uniform thickness and dewetted solutions in the absence of gas flow. We also predicted that for a relatively weak gas flow the liquid film may dewet in the centre but then heal and return to the uniform thickness state after the gas flow is switched off, as we indeed observe for q g = 0.525 slpm. For q g = 0.6 and q g = 0.7 slpm, we observe that the liquid film remains dewetted even after the gas flow is switched off, in agreement with the theoretical prediction of § 4.3.1. This is also confirmed using Gerris simulations.
For h 0 = 1 mm, the theoretical prediction was that the liquid film should heal after switching off the gas flow, no matter how strong the gas flow was. This is confirmed in figures 15, where we used the gas flow rates q g = 0.7, 0.8 and 0.9 slpm (see the black, red, and green lines, respectively). We switched off the gas flow at t = 0.2 s. The solid parts of the curves correspond to the gas flow switched on and the dashed parts correspond to the gas flow switched off. In all the cases, the liquid film returns to the uniform thickness state. Finally, we analyse dewetting of a water film of thickness h 0 = 1 mm for the gas flow rate q g = 1 slpm but for different equilibrium contact angles. We consider θ eq = 15 • , 30 • , 45 • , 60 • , 90 • , 135 • and 175 • in figure 16, illustrating interface profiles at time t = 3 s, at which steady states are reached. These results were obtained using Gerris, as for contact angles greater than 90 • our COMSOL implementation and thin-film model are not suitable. Videos showing the evolution of the water film for θ eq = 15 • and 135 • are included in the supplementary material. It can be observed that the approach to a steady state is non-monotonic in both cases, described instead by an oscillatory behaviour. In figure 16, we can observe that for larger contact angles the radius of the dry spot is larger, as expected. For θ eq = 15 • , we have confirmed the dry spot in the centre heals after the gas flow is switched off (as for θ eq = 30 • as we discussed above). However, for the larger equilibrium contact angles the liquid film remains dewetted when the gas flow is switched off.

Conclusions
We have examined both experimentally and theoretically the flow arising from a gas jet (air) impinging axisymmetrically on a liquid (water) in a cylindrical beaker. The calculations were carried out using two models that involved DNS (COMSOL and Gerris) and a third based on a thin-film approximation. The deformation of the liquid surface, determined experimentally when a steady state had been reached, was shown to be in good agreement with the DNS results for all water depths, when the gas flow rate was low, and in the case of a thin liquid film with all three models. It was also shown that interface shapes and streamline patterns calculated using the thin film model were in agreement with DNS for film thicknesses much larger than those for which the thin-film approximation was strictly valid.
Experiments were used to determine interface shapes both in the steady state and during the development of a steady flow after the gas jet was switched on. The contact angle between the liquid and its container was also measured experimentally and this was used as an input into the theoretical models. In the thin-film case, when the gas flow rate was high enough, dewetting of the film from the surface occurred. Using the experimentally measured contact angle of 30 • and a precursor thickness of h eq = 0.1 mm, a parameter determined from the disjoining pressure incorporating long-range and shortrange intermolecular forces and used in the COMSOL formulation and the thin-film equation, the conditions for dewetting determined by all three models were found to be in good agreement with experiment.
Dewetting was also investigated using linear stability analysis of various steady-state solutions of the thin-film model for a range of values of initial film thicknesses, contact angles and values of h eq . This analysis identified the various regimes in which dewetting could occur in agreement with the DNS and thin-film models. Regimes where the liquid would remain in its dewetted state or heal after the gas jet was turned off were identified.
For thicker films, the agreement between the models was less good for the time dependent flow before the final steady state was achieved, although the agreement for the steady states was still good. DNS results feature decaying interfacial oscillations/waves which were not present in the thin-film model, possibly due to the neglect of the inertial terms in the thin-film approximation. Experiments also showed the oscillations persisting for longer times than those predicted by the models (a video of an experiment for a 5 mm-thick water film and a gas jet of flow rate 1 slpm showing interfacial oscillations is included in the supplementary material). Indeed, preliminary experiments and DNS carried out over a range of gas flow rates for thicker films show that in some cases the oscillations do not decay at all, resulting in 'self-sustained oscillations' instead. The conditions for this to occur will be the subject of a future study.
Finally it should be pointed out that in the context of falling liquid films, accurate reduced-order models taking into account inertia have been developed using, for example, the weighted integral-boundary-layer approach (see e.g. Ruyer-Quil & Manneville 2000; Kalliadasis et al. 2011;Tseluiko & Kalliadasis 2011;Denner et al. 2016). Derivation and analysis of such models in the present context is left as a topic for future investigation.
We acknowledge the support of the UK Fluids Network (EPSRC Grant EP/N032861/1) in funding participation in a Special Interest Group meeting and a week-long Short Research Visit of RC to Loughborough University that shaped the direction of this work. CJO would like to acknowledge the Adventure Mini-CDT on Gas-Plasma Interactions with Organic Liquids at Loughborough University for the PhD studentship.