Linear and nonlinear instability in vertical counter-current laminar gas-liquid flows

We consider the genesis and dynamics of interfacial instability in gas-liquid flows, using as a model the two-dimensional channel flow of a thin falling film sheared by counter-current gas. The methodology is linear stability theory (Orr-Sommerfeld analysis) together with direct numerical simulation of the two-phase flow in the case of nonlinear disturbances. We investigate the influence of three main flow parameters (density contrast between liquid and gas, film thickness, pressure drop applied to drive the gas stream) on the interfacial dynamics. Energy budget analyses based on the Orr-Sommerfeld theory reveal various coexisting unstable modes (interfacial, shear, internal) in the case of high density contrasts, which results in mode coalescence and mode competition, but only one dynamically relevant unstable internal mode for low density contrast. The same linear stability approach provides a quantitative prediction for the onset of (partial) liquid flow reversal in terms of the gas and liquid flow rates. A study of absolute and convective instability for low density contrast shows that the system is absolutely unstable for all but two narrow regions of the investigated parameter space. Direct numerical simulations of the same system (low density contrast) show that linear theory holds up remarkably well upon the onset of large-amplitude waves as well as the existence of weakly nonlinear waves. In comparison, for high density contrasts corresponding more closely to an air-water-type system, although the linear stability theory is successful at determining the most-dominant features in the interfacial wave dynamics at early-to-intermediate times, the short waves selected by the linear theory undergo secondary instability and the wave train is no longer regular but rather exhibits chaotic dynamics and eventually, wave overturning.

(Dated: July 17, 2015) We consider the genesis and dynamics of interfacial instability in gas-liquid flows, using as a model the two-dimensional channel flow of a thin falling film sheared by counter-current gas. The methodology is linear stability theory (Orr-Sommerfeld analysis) together with direct numerical simulation of the two-phase flow in the case of nonlinear disturbances. We investigate the influence of three main flow parameters (density contrast between liquid and gas, film thickness, pressure drop applied to drive the gas stream) on the interfacial dynamics. Energy budget analyses based on the Orr-Sommerfeld theory reveal various coexisting unstable modes (interfacial, shear, internal) in the case of high density contrasts, which results in mode coalescence and mode competition, but only one dynamically relevant unstable internal mode for low density contrast. The same linear stability approach provides a quantitative prediction for the onset of (partial) liquid flow reversal in terms of the gas and liquid flow rates. A study of absolute and convective instability for low density contrast shows that the system is absolutely unstable for all but two narrow regions of the investigated parameter space. Direct numerical simulations of the same system (low density contrast) show that linear theory holds up remarkably well upon the onset of large-amplitude waves as well as the existence of weakly nonlinear waves. In comparison, for high density contrasts corresponding more closely to an air-watertype system, although the linear stability theory is successful at determining the most-dominant features in the interfacial wave dynamics at early-to-intermediate times, the short waves selected by the linear theory undergo secondary instability and the wave train is no longer regular but rather exhibits chaotic dynamics and eventually, wave overturning.

I. INTRODUCTION
Vertical two-phase flows of thin liquid films sheared by a counter-current gas are prototypical for many technical applications, such as absorption and distillation (using structured packing), evaporation and condensation. In these applications, not only mass and heat transfer but also operational limits are closely linked to the prevailing hydrodynamics. The flow in the two phases, in turn, is largely determined by the interactions between gas and liquid at its interface. Although gas-sheared liquid films have been part of active research on fundamental and practical level for several decades, the rich interfacial dynamics are still not fully understood. These rich dynamics emerge as the liquid interface becomes unstable, leading to the development of waves. Depending on the flow rates, these waves can form highly complex structures, which may give rise to wave break-up, ligament formation and droplet entrainment.
One of the first to investigate a channel flow with two superimposed fluid layers from a theoretical point was Yih [1], who used asymptotic expansion to solve the Orr-Sommerfeld eigenvalue problem associated with the temporal linear stability in the long wavelength limit for equal densities and layer thicknesses. He found that viscosity stratification alone can cause interfacial instability at arbitrarily small Reynolds numbers, which is also referred to as Yih mechanism. Yiantsios [2] extended the linear stability analysis by accounting for short waves as well as effects due to surface tension and gravity. They observed that, depending on the choice of parameters, the flow is receptive to a short-wave instability at low Reynolds numbers and, moreover, to a shear mode instability (Tollmien-Schlichting mechanism) for sufficiently large Reynolds numbers. To classify the various types of instabilities arising in parallel two-phase flow, Boomkamp et al. [3] analysed by which mechanism energy is transferred from the primary flow to growing disturbances, thereby verifying that both the Yih and shear mode are important routes to interfacial instability. In fact, a combination of these two mechanisms, also referred to as internal mode, represents a further possible route. An established method to obtain the full spectrum of eigenvalues in the Orr-Sommerfeld problem is the Chebyshev collocation method [4]. This method can also be used to determine whether the interface is convectively or absolutely unstable [5,6], a characteristic that is of importance for the operation of processes utilizing shear flows. An extensive review on the analysis of spatially developing flows can be found in Reference [7].
Linear stability analysis is an effective and proven technique to gain valuable insight in the genesis of interfacial instability. However, its validity is limited to disturbances with infinitesimal amplitude. As these disturbances grow, nonlinear effects become important and have to be taken into account. Due to the complexity of nonlinear stability only few general theories exist [8]. Nonetheless, a variety of modelling techniques have been proposed to describe the development of the interfacial waves up to finite amplitude. Some of these techniques impose a priori assumptions about the wave dynamics under investigation, like long-wave or lubrication approximation, resulting in model equations such as the Kuramoto-Sivashinsky equation or depth-averaged integral equations. A large number of studies on the nonlinear dynamics of interfacial flows are based on these kind of models [9][10][11][12]. However, given that their range of applicability is generally not known in advance, they may at times be prone to ambiguous results [13]. This is especially the case for fairly "punishing" flow regimes, involving large pressure fluctuations and potentially large-amplitude waves, for which there is a major necessity to gain fundamental understanding. On the other hand, weakly nonlinear theories based on either the Stuart-Landau or the Ginzburg-Landau equations dispense with assumptions that cannot be confirmed a priori and are capable of matching experimental observations [14][15][16]. Such "pure" weakly nonlinear theories are therefore appropriate to complement direct numerical simulations (DNS) of the full Navier-Stokes equations, which, in turn, are guided by basic Orr-Sommerfeld (OS) analysis, to study interfacial instability in a rigorous manner.
The aforementioned techniques have helped to shed light on the genesis and development of interfacial instability in shear flows and a good understanding of the mechanisms at play has been gained. However, a substantial amount of the available literature is dedicated to horizontal flows or flows down an inclined plane. Flow dynamics specific to a vertical configuration have not received as much attention even though the same methods are applicable. Phenomena related to this configuration, such as partial liquid flow reversal due to counter-current gas (flooding), influence the design, optimization as well as operation of technical processes to a great extent and have hitherto mainly been investigated experimentally. We therefore try to help elucidate the rich dynamics of vertical counter-current gas-liquid flows from a more theoretical standpoint using the (semi-)analytical methods and simulation techniques described above.
In this work we consider the laminar-laminar case only, although the ideas and results contained herein can be extended to turbulent gas streams. The most accurate methodology appropriate for such an enhanced level of complexity is full-scale DNS, which is the target of future work. However, short of full-scale DNS, a quasi-laminar model may be assumed for the linear stability of the two-phase flow [6,17] or, alternatively, a weighted residual integral boundary layer method (WRIBL) model [18][19][20], which also models nonlinear interfacial waves. These approaches both have their own advantages and shortcomings, and the aim of future work will be to confirm these reduced-dimensional modelling approaches with evidence from accurate DNS, towards which the present work is an initial contribution.
This work is organized as follows. We give a description of the investigated problem and outline the employed methods in § II. Results regarding temporal stability of the system are discussed in § III. Spatio-temporal behaviour with respect to absolute/convective instability of the linearized dynamics is presented in § IV. Section V discusses nonlinear wave dynamics. Finally, we give concluding remarks in § VI.

II. PROBLEM DESCRIPTION AND COMPUTATIONAL METHODS
In this analysis we consider the dynamics of a gas-liquid flow in a vertical channel, described schematically in Fig. 1. The two continuous phases are separated by an, initially, flat interface. A pressure gradient dp/dx > 0 in vertical direction counteracts gravity. We investigate cases in which the balance between gravity and pressure gradient gives rise to counter-current flow, with gas flowing in the direction of decreasing pressure and liquid flowing in the direction of gravity. Figure 1 also shows the development of a linear small-amplitude wave at the interface. Typically, the evolution of interfacial waves depends sensitively on the details of the mean flow.
Both fluid layers exhibit steady, spatially uniform, laminar and incompressible flow along the vertical rectangular channel. To describe this two-dimensional flow, we use a Cartesian coordinate system, (x, z), in which the flat interface is located at z = 0 and the confining channel walls are located at z = −d L and at z = d G , respectively. Within these boundaries, the fully developed liquid and gas layer occupy the regions −d L ≤ z ≤ 0 and 0 ≤ z ≤ d G , respectively.

A. Base flow and linear stability analysis
With the above mentioned conditions, the Navier-Stokes equations describing the fluid flow in both phases reduce to standard balances between pressure, viscous and gravitational forces, which are subject to no-slip condition at the channel walls, z = −d L and z = d G , as well as continuity of tangential stress at the interface, z = 0. To write the governing equations in nondimensional form we introduce the following dimensionless variables (without tildes) and scalings: where x = (x, z) and u = (u, w) are the coordinate and velocity vector, H is the channel height, V p is an inertial pressure scale, dp/dx is a positive pressure gradient, τ int is the interfacial shear stress, V * ,int is the gas side interfacial friction velocity and t denotes time. Further, the following dimensionless parameters arise: Here, µ j is the dynamic viscosity and ρ j the density of the respective phase (j = L, G), whereas δ j is the relative thickness of the respective fluid layers. The Reynolds numbers Re p , Re g and Re τ , in turn, relate to the applied pressure drop, to gravity and to interfacial shear. A Weber number W e accounts for surface tension. With this rescaling, the velocity profile for the undisturbed base flow in nondimensional form reads The ratio of pressure and gravity Reynolds number Re p /Re g results in a Froude number, which, in combination with the pressure scale in Eq. (1), represents a measure for the effect of applied pressure drop relative to gravity acting on the gas layer: To gain insight into the development of small disturbances η(x, t) centred around the flat interface z = 0, we examine the linear stability of the interface in terms of a standard Orr-Sommerfeld-type analysis (a comprehensive formulation of this analysis can be found in Appendix § 2. In this approach, the governing equations are linearized around the base flow, Eq. (3), and infinitesimally small perturbations in the associated streamfunction φ (z) are assumed to have wave-like solutions of the form Φ(x, z, t) = e i(αx−ωt) φ(z), where α = α r + iα i is the complex wavenumber and ω = ω r + iω i is the complex angular frequency (both dimensionless understood). The imaginary part of these quantities, α i and ω i , denote the spatial and temporal growth rates, respectively. If ω i > 0 the interface is considered temporally unstable and developing (initially) sinusoidal waves propagate with the phase velocity v p = ω r /α r .
We solve the underlying generalized complex eigenvalue problem numerically by employing a standard Chebyshev collocation method [4] (full description given in Appendix § 3), thereby adjusting the number of collocation points until convergence is reached.

B. Nonlinear direct numerical simulation
To capture and analyse the development of the flow beyond the linear regime, we use the two-phase flow solver presented in Reference [21] for direct numerical simulation of the full Navier-Stokes equations. This in-house solver is level set method based and uses the continuum surface force (CSF) formulation to model surface tension effects [22,23]. In level set formalism, the governing equations read as Here, φ (x, t) is the level set function indicating the phase in which a point x lies (liquid phase for φ < 0, gas phase for φ > 0). Hence, the zero level set, φ (x, t) = 0, represents the interface η (x, t). The level set function also determines the unit vectorn normal to the interface and the interface curvature κ in Eq. (5d). As the level set function differentiates between the two phases, it is also used to identify the respective density and viscosity through the expressions ρ = H (φ) + r (1 − H (φ)) and µ = H (φ) + m (1 − H (φ)). The function H (φ) is a regularised Heaviside function, which is smooth in an interval [− , ] around the interface with set equal to 1.5 times the grid spacing. This interval also supports the regularised delta function δ (φ) = dH (φ) /dφ. To discretize Eq. (5), we use an isotropic MAC grid with a spacing that resolves the undisturbed liquid film with at least 30 grid points. Additionally, all simulations have been checked for convergence. On the implemented grid, vector quantities are defined at cell faces and scalar quantities are defined at the respective cell centres. A third-order Adams-Bashforth scheme is used to treat the convective derivative, while the momentum fluxes are treated in a flux-conservative fashion employing a (semi-implicit) combination of Crank-Nicholson and third-order Adams-Bashforth method [24]. Pressure and the associated incompressibility of the flow are treated using standard projection method [25]. We use a combination of Jacobi's method and successive over-relaxation on a red-black scheme to evaluate the predictor and corrector step. Moreover, a third-order (fifth-order accurate) WENO scheme [26] is used to advect the level set function φ, which is subsequently reinitialised applying a Hamilton-Jacobi equation and the algorithm formulated in Reference [27].
At the domain boundaries we apply standard no-slip and no-penetration conditions at the confining walls, z = 0 and z = H (note that the coordinate system underlying the Eqs. (5) is shifted by −δ L compared to the one shown in Fig. 1), as well as periodicity in x-direction. The pressure is decomposed as p =p + (dP/dL) x, wherep satisfies the periodic boundary conditions in x-direction and dP/dL is a positive, constant, dimensionless pressure gradient. Solving the standard force balance in both phases (Eq. (A.1) and (A.5)) numerically gives the initial velocity field (base flow).
The solver is implemented in Fortran 90 using MPI to decompose the computational domain ( Fig. 1) in x-direction. This parallelization scheme makes efficient use of the architecture of the UK's "Advanced Research Computing High End Resource" (ARCHER, http://www.archer.ac.uk) on which we run our high resolution simulations.

III. TEMPORAL STABILITY ANALYSIS
In this section we restrict ourselves to the study of the temporal evolution of an (initially) infinitesimally small perturbation of the liquid interface. This temporal framework provides deep insight into the onset of interfacial waves and, thus, the complex dynamics of vertical films sheared by counter-current gas flows. It further enables us to map flow regimes typical for such systems.
We analyse the temporal stability for two distinct cases: • High density contrast: we assume a liquid density of ρ L = 1000 kg m −3 , corresponding to a gas-liquid flow typified by an air and water combination. We demonstrate below that the large density contrast leads to a complicated characterization of the instability, consisting of competing and coalescing linearly unstable modes.
• Low density contrast: we assume a liquid density of ρ L = 10 kg m −3 . The low-densitycontrast scenario is popular in the simulation literature [28][29][30] for a number of reasons. In the present context, it corresponds to a system without mode competition. This provides a "clean" database of linear stability results, which can be used as an unambiguous benchmark for direct numerical simulations.
For both considered cases, we further assume a liquid dynamic viscosity of µ L = 500 · 10 −6 Pa s and a surface tension of γ = 1 · 10 −3 N m −1 , while the dynamic viscosity and density on the gas side are µ G = 10 · 10 −6 Pa s and ρ G = 1 kg m −3 , respectively. The flow is confined by a channel of the height H = 0.01 m. These values result in density ratios of r = 1000 and 10, a viscosity ratio of m = 50 as well as a gravity Reynolds number of Re g = 313 and leaves the relative film thickness δ L and the Froude number F r, which can be related to liquid and gas flow rates, respectively, as the remaining parameters to determine the two-phase flow. Throughout the linear stability analysis, i.e. both density ratio cases, we consider δ L ∈ [0.02, 0.14], whereas the Froude number is varied within the interval F r ∈ [1.05, 13] for the case of high density contrast. This corresponds to an absolute film thickness d L ranging from 0.2 · 10 −3 m to 1.4 · 10 −3 m and an applied pressure drop dp/dx in the range of 10.8 Pa m −1 to 1657.9 Pa m −1 . Within this parameter space, we determine the temporal dispersion relation ω temp = ω (α r , α i = 0) numerically on a grid with ∆δ L = 0.005 and ∆F r = 0.25. The associated eigenvalue problem, Eq. (A.19), is solved for α r ∈ [0.05, α c ] with ∆α r = 0.05, where α c is the cut-off wavenumber beyond which ω temp i < 0. In the low-density-contrast case, we apply the interval F r ∈ [1.05, 1.55], corresponding to dp/dx ∈ [10.8, 23.6] Pa m −1 , with ∆F r = 0.025 and determine the dispersion relation for α r ∈ [0.05, 15] with ∆α r = 0.01. For each parameter set (δ L , F r) we further identify the pair α temp m , ω temp i that maximizes ω i as the linearly most unstable mode.
For both density ratio cases, our analysis shows ( Fig. 2) that the temporal growth rate of that mode always attains positive values, which means the interface is inherently unstable within the investigated parameter spaces. However, the trends for the growth rate of the with increasing Froude number in the case with high density ratio, a valley stretching out parallel to the horizontal axis at a Froude number of around 2.25 becomes evident in the δ L -F r-plane (Fig. 2a). Moreover, a plateau-like region of nearly equal temporal growth rates appears at about F r = 6.00 spanning all considered film thicknesses. Further insight into the nature of the system can be drawn from the length of the most unstable wave, λ m = 2π/α temp m , and its extent relative to the respective film thickness, λ m /δ L . Figure 3a shows the change of the latter quantity in the parameter space for density ratio r = 1000, revealing the complex interfacial dynamics of the flow. Although the interface is susceptible to a wide range of linearly most unstable waves, short-wave instability is predominant in this case, which is a consequence of the low surface tension considered herein [2]. The rich dynamics observed for the high-density-ratio case can be explained by looking at the dispersion relation of representative points in the parameter space (Fig. 4). It becomes apparent that multiple unstable modes are active, which are coalescing and competing with each other. It is this mode interaction that causes the complex structures presented above. In contrast, the case with low density ratio shows mostly one unstable mode throughout the considered parameter space. Even though further unstable modes are observed, these modes are very weak and appear only occasionally. Hence, for this case the system behaviour is less intricate (Fig. 2b and 3b).
In order to identify the character of these unstable modes and hence the driving force of the instability in these scenarios, we characterize the energy transfer from the base flow by decomposing the disturbance kinetic energy into production and dissipation terms. However, as both presented parameter spaces are quite extensive, we restrict ourselves to the detailed analysis of the representative scenarios listed in Table I. The mode characteristics for all other parameter sets are deduced by inspection.

A. Energy budget
Following the approach in Reference [3], the rate of change of kinetic energy in both phases KIN L,G can be decomposed into where all terms are spatially averaged contributions of the respective fluid to the overall energy budget of the disturbance. The terms DISS L,G represent energy losses due to viscous dissipation, whereas REY L,G denotes wave-induced Reynolds stresses transferring energy between the base flow and the perturbation in the bulk of the two phases. Lastly, N OR and T AN describe the work done by normal and tangential stresses at the interface. The energy budgets for the investigated scenarios are given in Table II, where the individual terms have been scaled by the total rate of change of kinetic energy (KIN L + KIN G = 1).
Decomposition of the energy budget in scenario T1C (Fig. 4f) uncovers three unstable (active) modes: one interfacial mode as well as one shear mode in each fluid layer. The first mode (α = 40.05) originates from the viscosity and density contrast of the two fluids, which give rise to velocity and stress disturbances tangential to the liquid interface driving the instability ("viscosity-gravity-induced instability"). Across the entire parameter space of the high-density-ratio case, this interfacial mode is active (Fig. 4) and, generally, its maximum growth rate ω temp i as well as the associated wavenumber α temp m increase for increasing Froude number yet slightly decrease with increasing relative film thickness, see Fig. 2a and 3a.
The shear mode, or Tollmien-Schlichting mode, in the gas layer is also active for all values of δ L above a minimum Froude number ranging from 1.46 to 2.34 for increasing film thickness. Above this threshold, the mode becomes stronger with increasing gas flow, though thicker films exhibit slightly lower growth rates. The opposite trend is observed for the wavenumber at maximum growth rate, which shifts towards lower values with increasing F r and generally higher values for thicker films. Contrasting the behaviour of the shear mode in the gas layer, the liquid layer shear mode is only active for thick enough films (δ L > 0.067). While this threshold value increases with increasing Froude number, the maximum growth rate of this mode drops, thereby rendering it less unstable. Further complexity is added to the system dynamics by the fact that the energy contribution to this particular mode due to tangential stresses at the interface becomes larger for increasing F r and reaches significant values (T AN ≈ 1). Eventually, T AN becomes the dominant term but with REY L still being of importance to overcome the restoring effects in both phase. This combination of energy sources driving the instability is also known as internal mode.
Apart from this change in mode identity, there is a further internal mode emerging at Froude numbers ranging from 5.21 to 5.55 for increasing δ L , e.g. scenario T2C or Fig. 4g. As this mode grows stronger with higher gas flow rates, the wavenumber α temp m at its maximum growth rate shift towards lower values. Furthermore, the mode starts to coalesce with the dominant interfacial mode, which leads to the formation of a second "hump" in that dispersion relation (Fig. 4h). For even higher F r, the coalescing internal mode undergoes a change in identity itself and forms a second active interfacial mode ( Fig. 4k and 4l). This "splitting" of the dispersion relation induces a region of rapidly changing wavelength of the linearly most unstable mode at intermediate Froude numbers (Fig. 3a). Although mode coalescence occurs throughout the entire range of investigated film thicknesses, splitting of the dispersion relation is not observed for low values of δ L (Fig. 4j). It is further worth mentioning that additional modes can become unstable, which are, however, mostly temporarily active and thus of minor significance.
The multitude of active and coexisting modes leads not only to mode coalescence but also to a certain amount of mode competition. In the high-density-contrast case this mode competition mainly occurs at low Froude numbers, where the growth rate associated with the interfacial mode is still comparatively low. The shear modes at low and high values of δ L , on the other hand, exhibit stronger growth and therefore constitute the dominant mode ( Fig. 4a and 4c). As F r increases, the interfacial mode picks up strength and supersedes the shear modes as the dominating mode ( Fig. 4d and 4f), which results in a jump of the wavenumber α temp m associated with the maximum growth rate. This jump explains the regions of relatively long waves at low and high values of δ L in Fig. 3a.
As mentioned above, the low density-ratio-case exhibits mainly one unstable mode, which is due to the density and viscosity contrast of the two fluids (Table II). Although both differences account for energy transferred towards the disturbed flow, the contribution related to the viscosity contrast dominates. Hence, this mechanism is, in general, consistent with the so-called Yih mode [1]. It is further apparent that the relative fraction of kinetic energy associated with the liquid phase increases with increasing Froude number. This rise, together with an enhanced energy dissipation, can be linked to more agitation in the liquid film as we will show in section V. The amount of energy dissipated in the gas phase, on the other hand, seems to drop, which is counter-intuitive for an increased F r. Yet, dissipation in the gas does increase in absolute terms but at a slower rate as the total kinetic energy. That, in turn, leads to the seemingly decreasing rate of dissipation in the gas phase. The same effect can be seen for N OR and T AN .
Another result worth noting is the change in sign of REY G as F r increases, turning wave-induced Reynolds stresses from an additional dissipative energy "sink" to an energy "source". This, in turn, can be explained by an unstable Tollmien-Schlichting mode appearing in the gas stream (scenario dT7F) that delivers energy to the interfacial instability, thereby suggesting a transition to turbulence in the bulk of the gas phase. Therefore, this particular scenario can conceptually not be regarded as a "pure" Yih-type instability any more but tends towards an internal mode. This positive contribution of wave-induced Reynolds stresses to the instability occurs throughout the entire parameter space but is shifted to lower Froude numbers for thicker liquid films. However, it has to be emphasized that the tangential stresses doing work at the interface are the dominant driving force of the instability in all presented scenarios for the low-density-ratio case (Table II). Across the en-tire parameter space the maximum growth rate ω temp i of this interfacial mode increases with increasing δ L and F r (Fig. 2b), whereas the associated wavenumber α temp m increases with increasing Froude number but decreases for thicker films, which in agreement with [12]. Overall, the system is predominantly receptive to long-wave instability under low-densityratio conditions (Fig. 3b). into two regimes: one in which developing waves exhibit a positive phase velocity and another in which the phase velocity is negative. With respect to the chosen coordinate system (Fig. 1), these regions correspond to waves propagating downwards and upwards, respectively. The vanishing v p at the demarcation (dashed lines in Fig. 5) between these regimes relates therefore to a standing wave. By crossing this line with increasing F r, liquid adjacent to the interface begins to move upwards with the developing wave. This (partial) upward flow of the liquid phase is considered as flooding [31]. The curve v p (δ L , F r) = 0, which is referred to as loading curve in chemical engineering, indicates the operational limit of vertical counter-current gas-liquid flows. Thus, the region below the loading curve constitutes the counter-current regime. The flow map for the case with high density ratio (Fig. 5a) reflects the complex interfacial dynamics described above. In regions with mode competition and changing dominant mode, the phase velocity changes at a great rate due to the jump of the wavenumber α temp m associated with the linearly most unstable mode. For low Froude numbers, this leads to a region of flooding amidst the counter-current regime at the thin film end of the parameter space, while a region of comparatively slow waves occurs at the high film thickness end. In contrast, the phase velocity changes smoothly in the low-density-ratio case due to the absence of mode competition (Fig. 5b).

B. Flow regimes
To illustrate the effect of the linearly most unstable waves on the onset of flooding, Fig. 5 also shows the curve of zero interfacial velocity (dotted line) for the undisturbed base flow of the respective case, which can be seen as the analogue to the actual loading curve. In the high-density-ratio case the two curves (with exception of the "island" region) almost coincide, i.e. linear waves have no significant influence on the loading curve (Fig. 5a). This can be explained by looking at the expression for the phase velocity, which can be decomposed (trivially) as v p = u 0,int + u 1 (α), where u 0,int denotes the interfacial velocity of the undisturbed base flow. Yet, a good approximation for u 1 (α) can be found in the dispersion relation for free-surface capillary waves [17,32], such that u 1 (α) ≈ ω cap /α, where ω 2 cap = γα 3 / (ρ L + ρ G ) (dimensional understood). Even in the absence of this particular approximation, it can be argued that u 1 (α) should be symmetric in the densities (at least for the present vertical configuration), and should be a monotonically decreasing function of the densities, as generally the inertia tends to affect the phase speed in a reciprocal manner, i.e. the larger the inertia the smaller the phase speed. Thus, it can be seen that that v p → u 0,int for high density ratio. This also implies that interfacial waves do alter the flooding curve in the case with low density ratio. The counter-current regime is significantly reduced by the instability for large parts of the parameter space (Fig. 5b). However, below a threshold of δ L = 0.036 the onset of flooding is pushed towards higher Froude numbers.

IV. ABSOLUTE AND CONVECTIVE INSTABILITY OF THE LINEAR DYNAMICS
Complementing the temporal analysis of the interfacial instability, we study the response of the system to a pulse on the liquid interface. In particular, we are interested whether this initially localized disturbance spreads both up-and downstream and eventually contaminates the entire flow, rendering the system absolutely unstable. In case the disturbance propagates away from its source, the flow is considered convectively unstable.
To determine the nature of the instability in this context, the complex dispersion relation D (α, ω) = 0, which is obtained by solving the eigenvalue problem of Eq. (A.17) numerically for a range of complex wavenumbers α = α r + iα i , has to be evaluated against conditions essential for absolute instability as outlined in Reference [7]: (i) a positive imaginary part of the angular frequency ω i,0 := ω i (α S ) > 0 at a saddle point α S in the complex α-plane, where α S solves dω/dα = 0, forms the necessary condition; (ii) to satisfy the sufficient condition, spatial branches α ± (ω) that originate from opposite halves of the α-plane have to coalesce at the saddle point α S , forming a genuine pinch point (this coalescence corresponds to the formation of a branch point/cusp at ω i,0 in the complex ω-plane [33]). Meeting both conditions will result in growth of the disturbance at its source with the absolute growth rate ω i,0 .
Although the described procedure seems straightforward, inspection of the results of a spatio-temporal Orr-Sommerfeld (ST-OS) analysis requires great care in order to avoid misinterpretation due to the complexity of the dispersion relation that can arise from the multivalued nature of the eigenvalue problem as well as specifics of the investigated problem, such as applied boundary conditions or multiple unstable temporal modes (see Fig. 6). Hence, to conclude the character of the instability correctly, it is necessary to examine the global topography of frequency over wavenumber for each relevant set of parameters. As we are interested in identifying convective/absolute instability (C/A) transition, the described method is not practical given the large parameter space considered herein. Instead, we make use of the analytical connection between temporal and spatio-temporal frequency presented in Reference [34]. The authors derived an expression for ω i (α r , α i ) based on analytical continuation of purely temporal stability properties, specified on the real axis, into the complex α-plane. Accounting for up to second-order terms, the quadratic approximation (QA) of the imaginary part of the complex frequency reads Similarly, we approximate ω temp (α r ) in the vicinity of the temporally most unstable mode by second-order polynomials. With the conditions ∂ω i /∂α r = ∂ω i /∂α i = 0 applied to Eq. (7), we can readily estimate the position of the saddle point α S and the absolute growth rate ω i,0 for a given set of flow parameters, requiring results from temporal linear stability analysis only. Even though this approach seems to circumvent the pitfalls associated with the full spatio-temporal linear stability analysis outlined above, one has to keep in mind that the quadratic approximation is, strictly speaking, only valid inside a disc of convergence with radius R, where R is the minimum distance from the centre of the complex Taylor series (α r = α temp m , 0) to the nearest singularity of ω (α). Like in the purely temporal framework, we study several scenarios in more detail (Table III). Using spatio-temporal Orr-Sommerfeld analysis and quadratic approximation as outlined above, we determine the saddle point α S as well as the corresponding absolute growth rate ω i,0 for all listed scenarios and further perform direct numerical simulations for the scenarios ST1 and ST4. In the following, we want to compare and discuss the obtained results for scenario ST4 in more detail.
As shown earlier already, the global topography of ω i in the complex α-plane (Fig. 6) is rather complex. First, the dispersion relation contains two saddle points, of which both may be dynamically relevant. The confinement of the flow by the channel walls has, furthermore, created a discrete pole on the imaginary axis [34], (α r , α i ) = (0, 3.34), which has implications on the character of the saddle point closer to that particular singularity. Lastly, the multivalued nature of the dispersion relation becomes apparent by the branch cut just below the real axis. Although these features make the final characterisation of the instability more difficult, the saddle point at (α r , α i ) = (4.87, 0.97) clearly appears as a result of the coalescence of spatial branches emanating from opposite half-planes. It is therefore also pinch point and contributes to spatio-temporal growth at a rate of ω i,0 = 0.5221. The positive value of α i indicates that the disturbance travels upwards, which is confirmed by DNS (Fig. 7a). On the other hand, the saddle point at (α r , α i ) = (0.62, 0.65) is not a pinch point. In fact, the spatial branch α + (ω) (above the saddle point) is a closed curve, whereas the α − (ω) branch does not (entirely) originate from the negative half-plane. Hence, this saddle point is dynamically irrelevant regarding absolute instability. As mentioned before, the singularity closest to the position of the temporally most unstable mode (α r = α temp m , 0) determines the disc of convergence in which the quadratic approximation can be applied with confidence. In scenario ST4, the confinement pole at (α r , α i ) = (0, 3.34) limits the outermost radius R of this disc to about 5.68. Equation (7) is thus convergent across the relevant section of the complex α-plane depicted in Fig. 6. The approximated position of the pinching saddle point as well as the corresponding growth rate agree well with spatio-temporal OS analysis, see Table III.
We further compare the impulse response for the parameters of scenario ST4 captured by DNS against linear theory (ST-OS). In contrast to simulations within the purely temporal framework, a Gaussian pulse with a height of 1 · 10 −3 and a standard deviation of 1 · 10 −1 is implemented on the interface as an initial condition to study the spatio-temporal nature of the flow. Growth of this perturbation at its source then constitutes absolute instability. However, the developing disturbances will inevitably contaminate the pulse source due to the implemented streamwise periodic boundary conditions. To delay this contamination sufficiently, the channel length is set to L x = 20. Using the pulse norm , the space-time plot in Fig. 7a shows the temporal evolution of the perturbations along the streamwise coordinate. It can be seen that, starting from its initial position at x = 16, the pulse is convected upwards, which is in accordance with the results from linear theory mentioned above. This upward motion triggers a pressure disturbance moving ahead of the pulse. However, this "shock" decays sufficiently before it re-enters the computational domain, thus avoiding early contamination of the instability source (see left-hand side of Fig. 7a). Further information from this plot is extracted in Fig. 7b, where the growth of the instability at its source is given. After an initial transient period, the growth rate of the instability stabilises and follows the value predicted by ST-OS and QA closely, therefore validating the DNS results. Thereafter, a second pressure shock develops (t ≈ 11.5) and travels upwards to eventually contaminate the pulse source at around t = 14.0. The excellent agreement that has been established between linear theory (ST-OS) and quadratic approximation in ST4 can also be observed for the scenarios ST2 and ST3 (Table III). For scenario ST1 we were not able to determine the absolute growth rate of the system by means of OS analysis and saddle point method due to the multivalued nature of the associated eigenvalue problem. In fact, the dynamically relevant saddle point appears in the lower half-plane below the branch cut arising near the real axis. To analyse this saddle point a laborious reconstruction of the corresponding part of the Riemann surface from discrete eigenvalues would have to be carried out. However, we avoid this procedure by using the quadratic approximation, which indicates scenario ST1 to be convectively unstable. A further direct numerical simulation confirms this (not shown).
Given these results, we apply the QA to identify the spatio-temporal nature of the system throughout the wide parameter space considered in the low-density-ratio case. The calculated growth rates ω i,0 are shown in Fig. 8, where the dashed lines demarcate the C/A boundaries. It becomes apparent that the system is absolutely unstable in almost the entire domain with exception of two narrow bands at the low Froude number and low film thickness end of the parameter range, respectively. In view of these results and those of § III B, it appears that C/A transition and the onset of flooding are not closely related for the case of density ratio r = 10.
Regarding the accuracy of these results, it has to be mentioned that the relative error (ω QA − ω ST-OS ) / |ω ST-OS | between QA and linear theory increases with decreasing film thickness and increasing Froude number. At the same time, the quadratic approximation generally underestimates the "true" value of ω i.0 as obtained by ST-OS, which leads to a larger regime of absolute instability (towards thin films) than displayed in the above figure. Finally, the mode coalescence and mode competition encountered in the case with high density ratio (e.g. Figure 4) makes it difficult to identify the nature of the instability in a spatio-temporal framework with both the present methods. To study this manifestation of the instability, alternative techniques, such as linearized DNS [6] (wherein the linearized equations of motion are solved as a Cauchy problem, but still within the framework the Orr-Sommerfeld linear operators), might be more suitable. A detailed investigation of this particular regime, i.e. high density contrast, will be left to future work.

V. NONLINEAR WAVE DYNAMICS
As temporal linear stability analysis permits the analysis of infinitesimally small interface perturbations only, we carry out direct numerical simulations to study the evolution of these disturbances up to finite amplitudes for both low and high density contrasts. To that end, we use the in-house solver described in § II B with streamwise-periodic boundary conditions and initial interface location where A 0 is some small initial amplitude (herein A 0 = 1 · 10 −3 ), N is the number of linearly unstable modes initialized, α n = n (2π/L x ) n is a wavenumber in streamwise direction and ϕ n is a random phase. Even though periodic boundary conditions do not reflect the behaviour of a real system, it is appropriate to consider this setup as it allows for a rigorous comparison of DNS results with linear theory as well as unambiguous results of the Fourier transform taken of the interfacial wave at finite times [21].

A. Low density contrast
To allow for rigorous comparison with linear theory, we perform direct numerical simulations for each flow regime of the low-density-contrast case (lower part of Table I). For that purpose, we set the wavenumber α n in Eq. (8) equal to the linearly most unstable mode α temp m of the respective scenario. Figure 9a shows the L 2 -norm of the wall normal velocity perturbation w over time for scenario dT4C and it can be seen that the disturbance grows with the rate ω temp φ j (z) at the crest of the developing wave also yields excellent agreement (Fig. 9b). Also the other scenarios follow the theoretical predictions in terms of growth rate in an equally convincing fashion (Table I).
It can further be seen in Fig. 9a that the wave initially grows exponentially before nonlinear effects gain importance at about t = 8.0. Beyond that point, wave growth slows down and eventually saturates, leading to a steady nonlinear wave of constant amplitude travelling along the interface. Figure 10a depicts the early stage evolution of the interface up to saturation for the counter-current scenario. It can be appreciated that the wave developing on the interface moves downwards as is predicted by linear theory. More insight into the flow features of the counter-current flow scenario is given by Fig. 11a. Plotted in a fixed frame of reference, a large anticlockwise rotating recirculation zone, positioned on the "downwind" side of the wave, is revealed in the gas phase. This vortex is sandwiched between the main flow of the gas and a thin layer of gas along the entire length of the interface, which is dragged downwards by the liquid due to interfacial shear stress. The liquid phase itself does not exhibit any major disturbances.
According to linear theory, the loading point of the system is reached by increasing the Froude number to F r = 1.179 (scenario dT5L, see Fig. 5b). At this point, the phase velocity of the wave vanishes as gravitational and lift forces on the wave hump balance each other, resulting in a standing wave. The evolution of this wave in its early stage is depicted in Fig. 10b. It can be apprehended that the initial upward motion of the wave slows down after it is saturated and eventually comes to a complete halt around the position last indicated in the plot. As the system evolves further the saturated wave begins to move downwards with the bulk of the liquid (not shown here). This is attributed to nonlinear effects which lead to an anticipated deviation from the system behaviour predicted by linear theory. However, it is worth mentioning that the deviation is relatively small (v p = 0.0103) and, hence, linear theory gives a good indication of the development and behaviour of the liquid interface for the loading scenario (dT5L). Moreover, the increased Froude number leads to a decrease of the mean liquid streamwise velocity, which, in turn, reduces the extent of the thin layer of downward gas flow along the liquid interface. Thus, the gas side vortex developing in the wave trough moves closer to the interface. The enhanced flow velocity on the gas side, on the other hand, stretches this eddy along the entire trough region (Fig. 11b). On the liquid side, contrary to the counter-current scenario, a small anticlockwise rotating vortex occurs in the wave body, which indicates a negative (upwards) streamwise velocity of the liquid at the interface in the vicinity of the wave crest. However, the bulk of the liquid, which refers to the region of the film between channel wall and wave trough, remains largely unaffected. By increasing the Froude number beyond the loading point, as in scenario dT6F, lifting forces on the wave body due to the enhanced gas flow overcome the gravitational forces, causing interfacial waves to travel upwards (Fig. 10c). Similar to the previous scenarios, the developing wave saturates and forms a coherent structure. However, with higher Froude numbers the general trend of increasing growth rate and decreasing amplitude becomes apparent. Furthermore, the change of direction in which the interfacial wave propagates leads to a more agitated liquid film, especially in the wave body (Fig. 11c). Compared to scenario dT5L, an extended region of liquid near the wave crest experiences upward movement with the wave body and becomes recirculated in an enlarged eddy. The bulk of the liquid, on the other hand, essentially remains unaltered also in this scenario.
After having discussed the influence of the Froude number, i.e. gas flow rate, on the gas-liquid flow in a qualitative manner, we now want to turn our attention to a qualitative description of the nonlinear wave dynamics. To understand the genesis of these, we compute the spectra of the interface height for each time step. The spectra of the first three harmonics are shown in Fig. 12a. Matching the rate given by Orr-Sommerfeld theory, the first harmonic (α 1 = α temp m = 3.99) grows exponentially fast at the beginning. At the same time, the higher harmonics (α 2 = 7.98, α 3 = 11.97) are linearly stable, as predicted by the same theory. However, these harmonics eventually also undergo exponential growth, whereat the n th harmonic grows at a rate nω temp i (α 0 ). As time goes by, the growth rate of all harmonics decreases and the amplitudes saturate simultaneously. It is thus apparent that the dynamics of the higher harmonics are "slaved" to the first harmonic. This temporal development of the interfacial spectra is perfectly consistent with weakly nonlinear theory [14]. The framework of the same theory allows to model the temporal development of the first harmonic using the Stuart-Landau (SL) equation [35,36]. Including up to quintic-order terms, the absolute value of the finite amplitude follows where µ 1 is twice the temporal growth rate 2ω temp i and µ 2,3 is twice the real part of the Landau coefficients [8]. These coefficients have been fitted numerically for the different scenarios discussed above and are listed in Table IV together with the root mean square deviation (RMSD) and R 2 -value of the best fit, illustrating the excellent agreement between theory and direct numerical simulations. The negative value of µ 2 confirms the existence of the supercritical bifurcation observed in all scenarios. For scenario dT4C, the relatively large value of µ 3 points towards a longer transition phase with decreasing growth of the wave amplitude leading from the initial stage of exponential growth to saturation in this scenario. In contrast, the wave in case dT5L, dT6F and dT7F reaches its saturated state faster, which is reflected by the vanishing of µ 3 in Eq. (9). The comparatively low R 2 -value in scenario T7F can be explained by an overshoot of the amplitude in the transition phase due to transient effects before the wave settles at a stable equilibrium amplitude.

B. High density contrast
To gain insight in the nonlinear wave dynamics under high-density-contrast conditions, we perform direct numerical simulations with the system parameters corresponding to Fig. 4e, scenario T3C (Table I). This system exhibits two distinct linearly unstable modes, one long-wave shear mode and a short-wave interfacial mode, which are accounted for by α REY G = 1.549 and α T AN = 40.27 in Eq. (8), respectively. It is due to the short-wave nature of the interfacial mode combined with the high inertia of the liquid phase that the flow characteristics near the interface (Fig. 13) differ significantly from those observed in the lowdensity-contrast case. Unlike the scenarios presented in the previous section, no prominent features, such as recirculation zones, emerge within the wave train. However, a vortex layer with the familiar "cat's eye structure" develops at the demarcation between the bulk of the gas and a thin gas layer dragged downwards by the liquid due to interfacial shear stress. Thereby, the vortices are pinched between two consecutive high-pressure regions, which are forming on the downstream side of the short-wave crests. Further snapshots indicate that this vortex layer is unstable to secondary instability. Hence, the wave form in Figure 13 should not be regarded as a quasi-steady state.
Similar to the low-density-ratio case, the nonlinear dynamics of the interfacial mode appear to be consistent with Stuart-Landau theory, albeit these dynamics develop faster due to the higher growth rate of the first harmonic (α T AN = 40.27). Evidence is provided in Figure 14. Consistency with the Stuart-Landau theory is especially clear in the time- frequency domain. In particular, in Figure 14b, a power-spectral density is shown (T corresponds to the duration of the simulation). A well-defined global maximum is observed at at ω r = 239.8, corresponding to the frequency of the most-dangerous TAN mode at α = 40.27 in the spatial domain. Successive maxima at ω r = n(239.8) with n = 2, 3, · · · indicate that the "overtones" of the TAN mode are slaved to the TAN mode itself.
The maxima in the power-spectral density function, although well defined, are by no means sharp. The broadening of the maxima here is a sign that the wave train is not strictly periodic, and that a very large number frequencies is present in the dynamics. Crucially, the broadness of the lines is a function of the simulation time: for longer simulations, more and more frequencies come into play (albeit that the same maxima remain dominant throughout). This indicates the onset of chaos. Thus, the Stuart-Landau theory manifests itself as the leading-order approximation to a chaotic dynamics, albeit that the wave train in Figure 13 is subject to secondary instability.
Other second-order effects are key to understanding the secondary instability. The vortex layer in Fig. 13a is not steady but breaks down at later times (Fig. 13b). In particular, a longwave perturbation to the vortex layer (wavelength on the domain scale) is seen to coincide with the significant growth of the long-wave shear mode. Therefore, it would appear that a complicated secondary instability sets in, involving a destabilization of the vortex layer by perturbations that are fed by the long-wavelength linearly unstable (but subdominant) mode. At this stage, individual waves steepen even further, up to the point where wave overturning is imminent (Fig. 13b).
One may furthermore recall that the strong presence of a shear mode in the system is a signal of the supercritical transition to bulk turbulence in one of the phases. This is not directly relevant in the present scenario, as the secondary instability of the laminar wave train sets in before this transition. In either case, it is clear that secondary instability may inhibit the operation of the system in a quasi-steady laminar state at high density ratios.

VI. CONCLUSIONS
We have presented a comprehensive study on two-dimensional laminar flow of a vertical liquid film sheared by laminar counter-current gas in a confined channel. This study tries to further elucidate the nature of interfacial instability in such two-phase flows using several complementary methods, namely Orr-Sommerfeld analysis, energy budget analysis as well as high resolution direct numerical simulations. Two sample systems have been selected for investigation: one with a high density contrast (ρ L /ρ G = 1000) and a second with a low density contrast (ρ L /ρ G = 10). In both cases, the same viscosity contrast (µ L /µ G = 50) and comparatively low surface tension (γ = 1 · 10 −3 N m −1 ) was used. In our study we focussed on analysing the influence of liquid film thickness and applied pressure drop on the development of interfacial waves.
Temporal linear stability analysis reveals that the liquid interface is inherently unstable for both cases. In the system with high density ratio short-wave instability is predominant, whereas the low-density-contrast case tends to favour long-wave instability. Furthermore, the instability is governed by a multitude of coexisting unstable modes (interfacial mode, shear mode in both phases, internal mode) under high-density-ratio conditions, where the latter two modes indicate the onset of turbulence in the bulk of one of the phases. In contrast, the low-density-ratio system exhibits only one unstable mode, which is consistent with the Yih mode [1] and a manifestation of the viscosity contrast.
Additionally, we use the phase velocity of the linearly most unstable mode to identify two distinct flow patterns: counter-current flow and flooding. In the former regime, the phase velocity is positive, which relates to interfacial waves propagating downwards. By contrast, flooding is characterised by upward travelling waves, leading to partial flow reversal in the liquid film. During the transition from counter-current to flooding regime, standing waves evolve at the loading point. In the case with high density contrast, we found a region of flooding amidst the counter-current regime, which is the result of mode competition (shear mode in the gas layer vs. interfacial mode).
We have further determined the nature of the instability in the spatio-temporal framework under low-density-ratio conditions. Using this information, we established a second flow map indicating the transition from convective to absolute instability (C/A). Besides standard Orr-Sommerfeld analysis, we also adopted the analytical connection between temporal and spatio-temporal growth rates in the linear regime as presented in Reference [34]. This approach, which is based on analytical continuation, circumvents possible difficulties in identifying the absolute growth rate that are associated with the multivalued nature of the eigenvalue problem and specifics of the problem at hand. Compared with OS analysis and DNS, this method shows good agreement and is therefore appropriate to accurately estimate the absolute growth rate of the instability. We find that the system is absolutely unstable in most parts of the parameter space considered herein with the exception of two narrow bands at low applied pressure drop and low film thickness, respectively. In the high-densityratio case, mode competition and mode coalescence between the multiple unstable modes hindered the mapping of convective and absolute instability in the respective parameter space. For such systems, linearized DNS [6] may allow for more conclusive results.
To assess the development of interfacial waves up to finite amplitudes, we perform direct numerical simulations for parameters within the established flow regimes of the low-densityratio case, using a level set method based solver that has been developed in-house. These simulations show excellent agreement with linear theory during the stage of exponential wave growth and confirm the determined flow patterns. The simulations also show saturation of the waves once nonlinear effects become important. A Fourier analysis reveals that the growth of the higher harmonics of the interfacial waves is coupled to that of the fundamental in a fashion which is consistent with weakly nonlinear theory. The growth of the first harmonic also agrees well with the Stuart-Landau model, thus underpinning the weakly nonlinear nature of the instability and, moreover, suggesting the existence of a supercritical bifurcation. Regarding high-density-contrast conditions, the dynamics of the (short-wave) interfacial mode appear to be similar to those observed under low density ratio. However, direct numerical simulations suggest that the emergence of an additional (long-wave) shear mode triggers a secondary instability, which leads to a chaotic wave train showing signs of imminent wave overturning. In summary, the combination of generic complementary (semi-)analytical and numerical methods presented herein yields a comprehensive and convincing characterization of interfacial instability in vertical counter-current gas-liquid flows. We are therefore confident that this rigorous methodology can be employed to further elucidate the dynamics of parallel shear flows in a wide range of technically relevant parameter regimes, such as flows with high viscosity contrast or non-steady gas flow. Beyond that, the outlined approach may be used as a starting point for the future study of heat and mass transfer phenomena in vertical gas-liquid flows. liquid film, this balance can be written as µ L d 2ũ 0 dz 2 − dp dx + ρ L g = 0, −d L ≤z ≤ 0, (A.1) whereũ 0 is the dimensional base flow velocity in the respective phase (tildes denote dimensional quantities). Equation (A.1) is subject to no-slip at the liquid side channel wall, z = −d L , and continuity of tangential stress at the interface,z = 0: Thus, integration of Eq. (A.1) yields as the velocity profile for the liquid film. The interfacial velocity, which constitutes one of the boundary conditions of the gas layer, reduces tõ The velocity profile for the laminar gas layer is derived analogous to the liquid layer. We therefore write the force balance as which is subject to continuity of velocity and shear stress at the interface: Applying these interfacial condition on Eq. (A.5) yields as the gas side velocity profile. To determine the interfacial shear τ int we apply the no-slip condition at the gas side channel wall,ũ 0 (z = d G ) = 0: In summary, the velocity profile of the undisturbed base flow in dimensional form reads (A.9) Applying the nondimensionalization scheme of Eq. (1) and (2) finally leads to the dimensionless form of the base flow velocity profile as presented in Eq. (3).

Numerical method
We solve this eigenvalue equation numerically by employing a standard Chebyshev collocation method [4] with an approximation for the streamfunction of the form where T n (·) is the n th Chebyshev polynomial of the first kind. We substitute this trial solution into Eq. (A.17) and evaluate the result at N L + N G − 6 Gauss-Lobatto collocation points. Together with the eight boundary and interfacial conditions, this yields N L + N G + 2 linear equations in as many unknowns. In matrix terms, we solve Lv = iωM v, (A. 19) where v = (a 0 , . . . , a N L , b 0 , . . . , b N G ) T . We use a standard linear algebra package (MATLAB R ) to solve this equation, thereby adjusting the number of collocation points (N L + 1, N G + 1) until convergence is reached.