Ultraefficient reduced model for countercurrent two-layer flows

We investigate the dynamics of two superposed layers with density contrast ﬂowing countercurrent inside a channel, when the lower layer is much thinner than the wavelength of interfacial waves. We apply a low-dimensional ﬁlm model to the bottom (heavier) layer and introduce a fast and efﬁcient method to predict the onset of ﬂow reversal in this phase. We study three vertical scenarios with different applied pressure gradients and compare the temporal growth rates of linear and weakly nonlinear waves to the Orr-Sommerfeld problem and to the weakly nonlinear theory, respectively. At the loading point, i.e., when a large wave hump stands at the interface, our spatiotemporal analysis shows that the system is absolutely unstable. We then present proﬁles of nonlinear saturated waves, pressure ﬁeld, and streamline distribution in agreement with direct numerical simulation. The reduced model presented here allows us to explore the effect of the upper-layer speed on the wave pattern, showing that the wave proﬁle is very sensitive when the mean ﬁlm thickness, rather than the liquid ﬂow rate, is maintained constant in the simulation. In addition, we show the strong effect of surface tension on both the maximum wave hump and the crest steepness before the loading point. Finally, we reveal how the nonlinear wave speed affects the vortex distribution within the lower layer by analyzing the stream function under different scenarios.


I. INTRODUCTION
Thin films sheared by a countercurrent flow have wide industrial applications, such as in absorption and distillation processes. During the past decades, this topic has been the subject of several experimental observations [1][2][3][4][5], aiming to elucidate the role of the flow pattern on the absorption rate. In addition, further to Yih's work on two superposed layers under viscosity contrast [6], several studies have investigated the linear stability of the resulting interfacial waves [7][8][9].
In the above-mentioned applications, the liquid film is often very thin compared to the wavelength of interfacial waves; this has encouraged the development of integral models in which the film dynamics is enslaved to thickness and flow rate only, with the main advantage of providing faster numerical simulations with respect to direct numerical simulation (DNS). Since the pioneering works of Benney [10] and Shkadov [11], these low-dimensional models have gained great accuracy. Particularly, through a weighted residual technique applied to the boundary-layer equations (WRIBL), Ruyer-Quil and Manneville [12] have addressed, for gravity-driven liquid films, the issues of blowup at moderate Reynolds numbers, typical of the Benney's model, and deviation from the long-wave stability threshold affecting Shkadov's model. Nowadays, low-dimensional models are also used for films under a constant shear [13] and for flows in a slowly varying duct [14]. Concerning the multiphase systems, the integral models have been first employed by Jurman and McCready [15] for liquid films sheared by a turbulent gas flow; later studies have extended the integral formulation also to the second phase, providing an evolution equation for the interface between two superposed layers [16,17]. However, these works follow those of Benney and Shkadov and thus experience the aforementioned shortcomings. Amaouche et al. [18] have instead applied the WRIBL model to both phases coflowing within a channel; nevertheless, they have accounted only for gravity-driven flows and their model cannot be used for countercurrent scenarios, as this work does. Another interesting approach has been developed by Tseluiko and Kalliadasis [19], where the WRIBL model is applied to the liquid phase for liquid films sheared by a countercurrent turbulent gas flow. However, from the gas side, the liquid film is modeled as a wavy wall and thus flows at comparable speed as treated in this work cannot be reproduced. Subsequently, their model has also been extended to coflowing channel flows [20]. Recently, Dietze and Ruyer-Quil [21] have developed a full two-phase WRIBL model with application to co-and countercurrent flows in narrow channels and then extended it to narrow tubes [22]. Nevertheless, using an integral model in the gas phase limits the thickness of the gas layer and its speed, resulting that large channels as considered in this work cannot be analyzed. Lately, Lavalle et al. [23] have adopted a first-order integral model to study laminar coflowing gas-liquid channel flows, which the present work is based on.
A peculiar feature of countercurrent flows is that under certain conditions they manifest flow reversal (or flooding), with dramatic consequences on the performances of technological devices. Various works have therefore aimed to predict the onset of flow reversal using experimental correlations (see review from [24] and later investigations [3,4,25]); in support, several studies have provided new insights and novel methodologies to analyze the flooding through linear stability [9], DNS [26], and low-dimensional models [9,19,21,22].
With this work, we provide an efficient and nonlinear method suitable for industrial uptake for predicting the onset of flow reversal, aiming to cut down computational cost and data with respect to DNS. We address the problem of countercurrent two-layer flows at low density contrast inside an inclined channel, assuming the lower layer thinner than the interfacial wavelength. The surface tension approaches the typical value of solvent extraction [27]. Our methodology originates from the work of Lavalle et al. [23], whose shallow water Arbitrary Lagrangian-Eulerian Navier-Stokes model (SWANS) consists of first-order depth-integrated equations applied to the lower layer and compressible Navier-Stokes equations for the upper layer. However, as it is, the SWANS model is unsuitable for countercurrent flows, because the shear stress exerted by the upper phase is negative and might cause a blowup of the Rusanov numerical flux when the characteristic velocity is negative. This flaw has been addressed in the current work by discretizing the term responsible for blowup with a central difference approximation and by changing the Rusanov flux accordingly. Furthermore, although it does not account for second-order streamwise viscous dissipation, this model provides good accordance with the Orr-Sommerfeld theory and DNS due to the smallness of the Weber number in the performed analysis. Finally, given that part of this study is based on the work of Lavalle et al. [23], we direct the reader towards it for further mathematical and numerical details.
The article is structured as follows: Section II introduces the low-dimensional model and the coupling technique with the upper phase. Section III shows the linear and nonlinear results and the comparison with the Orr-Sommerfeld problem, weakly nonlinear theory, and DNS. Section IV explores the influence of upper-layer speed and surface tension on the wave topology, as well as the role of the nonlinear wave speed on the vortex distribution. And finally, the conclusions and the perspectives of this work are discussed in Sec. V.

II. FILM MODELING AND COUPLING WITH THE UPPER PHASE
Consider a thin layer flowing under the gravity inside an inclined 2D channel, sheared by a countercurrent flow driven by an applied pressure drop exerting a shear stressτ i and a pressurep i at the interface (Fig. 1). The bottom layer flows due to gravity but also feels the opposite pressure gradient. The assumption that the thickness of the bottom layer is smaller than the wavelength 014001-2 of interfacial waves, i.e., ε =h/λ 1, yields that the variation in time and along the streamwise direction are smaller than those in the crossstream direction; it is thus found that∂ x ∂ t ∂ y , where the tilde defines dimensional quantities. Under these assumptions, the Navier-Stokes equations in the film turn into the boundary-layer equations. Their dimensionless forms up to O(ε), along with the corresponding boundary conditions, read The scaling is based on the uniform film thicknessh 0 and velocityŨ 0 , but other choices are possible. Consequently, Re =h 0Ũ0 ν −1 is the Reynolds number, Fr =Ũ 2 0 (gh 0 ) −1 the Froude number, and We = ρŨ 2 0h 0 γ −1 the Weber number. As a matter of fact, the use of the boundary-layer equations (1) as the beginning of the integral-model development allows us to relax the validity of such models [28]: it must be verified that ε Re 1 and Fr ≈ Re. Equation (1c) can be directly integrated by means of the condition (1e), leading to the hydrostatic pressure distribution where G(x) = −∂ x p i is the interfacial pressure gradient. After replacing (2) into the x-momentum equation (1b), we can perform the integration over the film thickness, which is the main ingredient of thin-film low-dimensional models. We obtain where q = h 0 u dy is the liquid flow rate, ∂ y u| 0 the wall shear stress, and = Re Fr −1 sin β + ReG(x). A closure model is required for the wall shear stress and the integral of squared velocity, and several ways exist to tackle this problem. We use an asymptotic expansion of the velocity field in terms of ε, where the leading order corresponds to the parallel flow and the first order to the O(ε) corrections. After rewriting the equations (1) at O (1) and O(ε) separately, we can find the velocity profile u (0) + u (1) , where the superscript (0) denotes the leading order while the superscript (1) the corresponding correction. The Appendix shows this procedure along with the resulting leading-and first-order velocity fields. The performed velocity profile is necessary to close the above-mentioned terms of (3b), keeping in mind that the right-hand side is one order greater than the left-hand side 014001-3 and requires O(ε) corrections. In other words, in order to fulfill the consistency, we use only the leading-order velocity (A1) to resolve the integral of squared velocity, while we need also the first-order velocity correction (A5) for the wall shear stress term. Details of this approach are given in Ref. [23]. Finally, the x-momentum equation (3b) becomes For numerical reasons, we gather all ∂ x h terms at the left-hand side inside P , while T contains the remaining correction terms of the wall shear stress; we also force the coefficient of q 2 /h equal to unity for the Galilean invariance [29], without losing any properties of consistency of the model. In this procedure indeed, we use the definition of the flow rate (A2). The terms P and T are defined as (derivations are given in Ref. [23]) The characteristic velocity e 2 = dP /dh = 2 Fr cos β must be positive to guarantee that the system of equations is hyperbolic (Toro [30]). While this is true for the majority of coflowing two-layer flows (see Fig. 4 in Ref. [23]), it is not always the case for countercurrent scenarios, as in this work. Specifically, the negative shear stress τ i and the resulting term 1/15 τ i h 4 in the expression (5) might cause e 2 to be negative. When this happens, the Rusanov flux f , computed at the generic cell face i + 1/2, manifests blowup because the characteristic velocity becomes complex. The Rusanov flux f reads indeed where v = (h,hU,hw) T is the conserved variable vector and f (v) = (hU,hU 2 + P ,hUw) T the corresponding flux, whereas w = h is the additional numerical variable that allows us to reduce the order of the discrete system [31]. Quantities v + i+1/2 and v − i+1/2 are the right and left states across the cell face in i + 1/2, respectively. To tackle the problem of blowup, we compute e 2 by using max[0,1/15 τ i h 4 ] at each time step. Note that this procedure is equivalent to discretize the term responsible for blowup with a centered scheme. With reference to the expression (7), the first part of the modified Rusanov flux remains then unchanged, while the second part is free of the complex contribution. By doing so, only the numerical viscosity of the scheme is changed (Toro [30]), while the discrete equations remain consistent with the physics described by the model (3a), (4)- (6). As it is, the modified model is tailored to countercurrent two-layer flows.
Nevertheless, in addition to the solution proposed above, two alternatives can be sought to get rid of the negative characteristic velocity e 2 : (i) the use of a fully implicit numerical scheme, which actually increases the complexity of the solver without guaranteeing a better structure than this one; (ii) the development of a three-equation film model, towards which the present work is meant to be extended in future, in line with the study of Richard et al. for gravity-driven falling films [32].
At this stage, the low-dimensional model (3a), (4)-(6) is coupled to compressible Navier-Stokes equations in the top layer. In dimensionless form those read  Layer where the film thicknessh 0 , the mean film velocityŨ 0 , and the uniform upper-phase densityρ 2 and pressurep 2 have been used as length, velocity, density, and pressure scales, respectively. The fluid is considered ideal; M =Ũ 0 /ã is the Mach number, Re 2 =Ũ 0h0 /ν 2 is the Reynolds number, while the Froude number Fr has the same definition as the bottom layer. In the limit of low speed flows, the viscous stress tensor reads T 2 2D 2 , where D 2 is the strain tensor. In the SWANS model, the coupling between Eqs. (8) and (3a), (4)-(6) assures the transfer of the interfacial stresses from the top to the bottom phase; vice versa is the transfer of the interfacial velocity and the shape of the interface. From a numerical point of view, this methodology makes use of an augmented system in the bottom layer accounting for an evolution equation of the surface tension to avoid scheme instabilities [31], as well as an accurate low-Mach scheme [33] with a moving mesh in the upper phase. The coupling technique is explicit and first order in time, while the boundary conditions are periodic. All the discretization techniques are fully developed in Ref. [23].

III. LINEAR AND NONLINEAR DYNAMICS OF THE COUPLED SYSTEM
Using the SWANS solver, we aim to study countercurrent two-layer vertical flows at low density contrast. Table I lists the physical properties of the two phases: density and viscosity ratios are chosen for numerical reasons, i.e., computational time of DNS, which we compare our model to. The performed direct numerical simulations are obtained by means of the TPLS flow solver, making use of the level-set method [8,34].
We run SWANS simulations varying the applied pressure gradient and the film thickness in accordance with the most unstable wavelength: we look at the sign of the displacement of the  [35]. The Froude number is now defined asFr 2 = p/λ/(ρ 2 g). traveling waves computed by the solver SWANS and cover the map in Fig. 2. This map quickly reveals the regions where the wave is moving upwards against the gravity. This phenomenon is crucial due to its potential towards leading to flooding, the operational limit of industrial grade absorption, and distillation columns. Furthermore, Fig. 2 has an immediate interpretation towards experiments and industrial procedures, because it links the applied pressure gradient driving the gas flow to the liquid flow rate, which can be recovered from the uniform film thicknessh 0 through the analysis of the base state (or from the leading order withh =h 0 ). In detail, three different situations appear in Fig. 2: a countercurrent saturated interfacial wave travels down with the gravity (DW for downwards wave), stands in a fixed position (SW corresponding to the so-called loading), or travels upwards driven by the top layer (UW as flooding). We have thus provided a nonlinear and highly efficient method to predict the onset of upward-traveling waves; a comparison with the standing-wave curve obtained through a solution of a corresponding Orr-Sommerfeld eigenvalue problem shows excellent agreement. However, the main advantage of our approach with respect to the Orr-Sommerfeld analysis is the straight extension to those configurations with nonflat walls and more general geometries. In addition, the low-dimensional formulation of the SWANS solver allows for a fast solution and returns much lower computational data in comparison to a full two-phase DNS simulation of the same kind. Among the SWANS simulations performed in Fig. 2, we focus on the three scenarios described in Table II. Noteworthy is that We is small and this explains why we use a first-order film model, given that the damping by viscous dissipation plays less because of the smoother interface. An initial perturbation of the uniform flow is applied in the form of h = h 0 [1 + sin(2πx/λ)] (and so for q), where the amplitude is = 10 −3 . For each test, Fig. 3 shows the temporal growth of the first three harmonics (with wave numbers k = 2π/λ, 2k, and 3k, respectively). We compare with good agreement the numerical growth of the first harmonic to the Orr-Sommerfeld problem. We also notice that the second and third harmonics are enhanced by nonlinear interactions governed by the fundamental harmonic; at this stage, they grow with rate proportional to the first harmonic [36], i.e., ω i{n} = nω i{1} , as also observed in our simulations. Remarkably, the model is able to capture the temporal growth of the harmonics, although their length is comparable with the film thickness, and again this must be due to the small We number.

014001-5
We further check whether the SW mode is connected to an absolute instability. Indeed, Vellingiri et al. [9] have investigated the instability of gas-liquid film flows showing that the generation of the standing wave occurs at the upper limit of the absolute instability region. Lately, Schmidt et al. [35] have demonstrated that at low density contrast the loading curve is always associated to an absolute instability. In order to inspect this tendency in the SW case, we study the response of the system subject to a localized disturbance by applying at the equilibrium a Gaussian pulse of the form f (x,0) = 2.5×10 −3 exp[−(x − c 0 )/2w 2 ], centered in c 0 = 0.7L with width w = L/40. However, the simulation with periodic boundary conditions requires a sufficiently long domain to avoid wave contamination; the length of the domain has been therefore chosen as ten times the wavelength associated to the SW scenario (see Table II). Figure 4 shows that after a sufficiently long time the pulse develops in both directions and the instability grows at the source, thus confirming that the system is absolutely unstable.
Moving back to the initial sinusoidal perturbation, the main wave saturates once the nonlinear effects begin to play in the system. The last two columns of Table II compare the speed c of the saturated waves between our model and the DNS; the small discrepancy is to be expected and has been also discussed by Ruyer-Quil et al. [37]. Figure 5 depicts instead the interfacial wave at different times for these three scenarios. As the waves saturate, agreement of the wave profile is shown between SWANS and DNS. The mismatch can be due to the omission of second-order elongation stresses in the film model and to the coupling strategy, which is quasistatic in the sense that the coupling is explicit in time. Figure 5 also shows that the amplitude of saturated waves slightly decays when increasing the applied pressure gradient [see Figs. 5(a)-5(c)]. This unusual behavior is explained because the applied pressure gradient has been mostly increased by reducing the imposed wavelength (i.e., growth of the frequency in a space-dependent experience), driving to a decrease of the amplitude, as observed by Liu and Gollub [38] for falling liquid films. Figure 6 represents the pressure field and the wall-fixed streamline pattern inside the lower layer and in the proximity of the interface. The streamlines are computed as a contour of the stream function = y 0 u dy, where the leading-order velocity (A1) has been used for the lower layer. We notice again agreement between the DNS [Figs. 6(d)-6(f)] and our model SWANS [Figs. 6(a)-6(c)], which captures well the different topology inside the heavier phase. We remark that the behavior in the lower layer extensively varies in the three flow scenarios; on the contrary, the upper-layer vortex generated by the wake behind the wave slightly changes because so does the wave profile among these three cases. Particularly, we observe a small tendency of our model to squeeze this vortex against the interface; again, this can be due to the coupling strategy, which is more sensible at the interface. Noteworthy is that the recirculation zones captured by our integral model in the wall-fixed reference frame are consistent with the long-wave approximation. Indeed, the shear stress τ i in the velocity profile (A1) might be negative due to the countercurrent flow, and as a consequence, the lower-layer velocity is not monotonic, as it happens for cocurrent film flows instead, permitting recirculation zones underneath main wave humps.

IV. ANALYSIS OF WAVE PATTERN AND VORTEX DISTRIBUTION
In this section, we describe the wave topology in the vicinity of the loading condition and explore the influence of both the upper-layer speed and the surface tension on the wave amplitude, speed, and steepness. Following Barthelemy et al. [39] and with reference to Fig. 7, we define the crest steepness as S c = aπ/λ c ; similarly, an analogous definition for the trough steepness S t is introduced. Finally, we investigate the velocity field within the lower layer with a particular focus on the genesis of recirculation zones occurring in the three scenarios previously discussed.

A. Influence of the upper-layer speed
Concerning turbulent gas-liquid flows, previous works [19,26] have shown that increasing the countercurrent gas speed leads to an increase in the amplitude of the wave when the Reynolds number of the film is fixed; meanwhile, the maximum and minimum of the wave profile increase and decrease in amplitude, respectively. With the present study, we supplement the analysis of the wave pattern by describing the influence of the upper-layer speed in the case of laminar two-layer flows at low density contrast. In addition, we show that the wave topology is very sensitive to the parameter that we fix in the simulation, e.g., equilibrium film thickness or mean liquid flow rate. Indeed, it is worthy to recall that fixing The gravity contribution of the flow rate, i.e., q g 1 = 1/3h 3 Re Fr −1 , corresponds to the flow rate of a vertical liquid film in a passive gas atmosphere. Therefore, fixing q g 1 in the simulation is the same as fixing the uniform film thickness.
In our work instead, the Reynolds number is defined through the effective liquid flow rateq 1 , namely, Re =q 1 /ν, whereby fixing Re is equivalent to keeping the liquid flow rate constant, which has a functional meaning for experimental reproduction. In what follows, we describe how the wave profile is sensitive to the upper-layer speed, starting from maintaining constant q g 1 and then q 1 . All the values of U 2 encountered in this section are meant as mean speed in the upper phase and are scaled with the mean speed of the SW case U SW = −2.69 m/s. Figure 8 shows that as the upper-layer speed increases, the saturated wave profile is strongly affected when q g 1 is fixed in the simulation (i.e., constant film thickness h 0 ). It is important to recall that the overall film flow rate changes with U 2 to accommodate a constant q g 1 . Figure 9(a) exhibits instead that the maximum (minimum) height slightly increases (decreases) when U 2 augments; meanwhile, the nonlinear wave velocity decreases [ Fig. 9(b)]. This behavior is in line with liquid-gas turbulent flows [19,26], although the rate of variation of the wave height is now smaller due to the different conditions adopted here: laminar flow, low density contrast, lower surface tension, and smaller wavelength. Also, increasing U 2 , the crest steepness S c increases whereas the trough steepness S t is slightly affected [see Fig. 9(c)].
On the other hand, the wave pattern changes extensively when the effective liquid flow rate q 1 is fixed while the upper-layer speed is increased. This corresponds to the case of a vertical film with fixed flow rate at different gas speeds and has a proper physical meaning. Figure 10 sketches the wave profile at different upper-layer speeds: unlike the previous case, we note that an increase of the countercurrent upper-layer speed barely modifies the shape of the saturated wave. Furthermore, Fig. 11(a) shows that the wave crest increases in amplitude when U 2 increases, and the minimum film thickness increases as well, unlike before. This happens because, at constant flow rate, the liquid decelerates when the countercurrent upper-layer speed increases, and the equilibrium film thickness 014001-10 increases accordingly, given that q 1 = hU 1 = const. Therefore, the increase of the maximum and minimum wave height is sustained by the increase of the main film thickness with U 2 . Conversely, the steepness of the crest and the trough behaves similarly as the previous case [ Fig. 11(c)]. Finally, the wave profile is more sensitive to U 2 when the equilibrium film thickness, rather than the liquid flow rate, is maintained constant.

B. Influence of the surface tension
In order to investigate the effect of the surface tension on the wave topology, the wavelength has been increased compared to the previous cases to avoid a fall in the stable region (surface tension stabilizes short waves). At p/L = 14.28 Pa/m, h 0 = 0.8 mm, and γ = 10 mN/m, we chose the wavelength of the most unstable mode corresponding to the loading curve, namely,λ = 45.9 mm. At this value ofλ, we explore the effect of increasing the upper-layer speed U 2 and the surface tension γ , maintaining q g 1 constant (i.e., film thickness h 0 constant). This choice is due to the previous analysis, which emphasized that the greatest variations of the wave pattern with U 2 occur when q g 1 is fixed rather than at constant flow rate q 1 . Noteworthy is that varying the surface tension must be consistent with the asymptotic velocity field (A5), which defines the dimensionless group = (h 0 /λ) 2 Re/We as distinctive for the change of the capillary term.
As the surface tension augments, Fig. 12 sketches the saturated wave profile at U 2 = 0.75, scaled with the mean upper-layer speed of the loading point U l = −4.19 m/s. Note that for this value of the upper-layer speed, the waves are moving downwards. We notice that an increase of the surface tension leads to a decrease of the amplitude of the saturated wave. This is also confirmed by Dietze [40] for vertical falling liquid films in a passive gas atmosphere. This trend allows us also to explain ]. In addition, increasing the surface tension γ leads to a decrease of both the maximum wave hump and the crest steepness, and this is more notable before the loading point (vertical bar in each plot). Indeed, before the loading point the effects of inertia and gravity overcome those of the countercurrent flow, and the liquid film is affected with γ as in a vertical falling scenario in a passive gas atmosphere [40]. Conversely, after the loading point only the nonlinear wave speed is affected by the surface tension, with a tendency of accelerating the wave downwards when γ increases. Finally, it is interesting to compare the results in Fig. 13 (solid lines only) with those in Fig. 9, where the wavelength is smaller. We find that increasingλ leads to an increase of the peak of the saturated wave before the loading point. This tendency is also in agreement with the vertical falling liquid films described by Dietze [40].

C. Influence of the nonlinear wave speed on the internal agitation
Lastly, we focus on the flow dynamics in the reference frame moving with the wave speed, where streamlines and trajectories coincide. This allows us to elucidate the regions of high transfer and is therefore relevant for industrial uptake. Figure 14 shows the streamlines in the reference frame moving with the wave speed for the DW [Figs. 14(a) and 14(b)] and UW [Figs. 14(c) and 14(d)] tests and compares with good agreement the model [Figs. 14(a) and 14(c)] to the DNS [Figs. 14(b) and 14(d)]. Note that since in the SW case the wave stands in a fixed position, the fixed and moving reference frames coincide. Remarkably, we notice two vortices in the DW scenario, and zero for the UW. In the SW case, instead we observe one single vortex that is clearly located at the wave hump [see also Fig. 6(b)]. Interestingly, the two vortices of the DW case are similar to those observed by Trifonov in countercurrent air-water flows [41]. In order to clarify the presence and the location of those vortices, we write the stream function in the lower layer, taking into account the wave speed c. Integrating the velocity profile (A1), one gets (the constant of integration vanishes for consistency with the cross-stream velocity v (1) ) Figure 15(a) sketches the stream function at the wave crest abscissa; a necessary condition for the presence of the vortices is that the stream function have a local minimum or maximum. We then realize that for c = 0 (SW case), only one extreme can be observed and is located next to the interface; if c < 0 (UW case), no vortices (or possibly a small one at the top) can be recovered; and when c > 0 (DW case), we can observe also a vortex next to the bottom wall, as confirmed by our simulations in Fig. 14. In Fig. 15(b), the stream function at the wave trough abscissa is sketched instead. In this case, the absence of local extrema suggests that no vortices exist in that region, although a small recirculation zone can be observed at the bottom when c < 0. Again, this is entirely confirmed by our simulations in Fig. 14, in accordance also with DNS. Interestingly, Eq. (10) highlights that the stream-function distribution between crest and trough changes thanks to the evolution along x of the shear stress, the pressure gradient, and the interface. For the sake of precision, note that since h(x) [and therefore the related quantities τ i (x) and G(x)] slightly changes among the three performed tests, the stream function of Fig. 15 corresponds to the SW case, to which we added a prescribed wave speed c. Summarizing, we reveal that the wave speed affects the agitation and the vortex distribution within the lower layer.

V. CONCLUSIONS
In conclusion, this work is based on the development and application of a first-order lowdimensional model to countercurrent two-layer flows in a channel. This model, coupled to compressible Navier-Stokes equations in the upper phase through a moving mesh technique, provides an efficient method to predict the onset of flooding. In presence of low density contrast and low surface tension, the temporal growth of linear and weakly nonlinear waves matches with the results of the Orr-Sommerfeld and the weakly nonlinear theory. An analysis of the response of the system to a localized Gaussian pulse shows that the loading point is connected to an absolute instability. Furthermore, the comparison of saturated wave profiles, pressure field, and streamline patterns shows agreement with the DNS. Using this coupled model, we then investigated the influence of upper-layer speed and surface tension on the wave pattern. We have found that increasing the upper-layer speed has a stronger effect when the liquid film thickness is maintained constant, rather than when the liquid flow rate is kept fixed in the simulation. When varying the surface tension instead, both the maximum wave hump and the crest steepness are affected before the loading point. Finally, we have given some insights regarding the vortex distribution within the lower layer, revealing that when the saturated wave is traveling downwards, more recirculation zones are to be expected.
The extension of this work would move towards the analysis of gas-liquid flows, taking into account the turbulence and then including mass transfer. Of high interest is also the issue of whether a low-dimensional integral approach is able to capture the mode competition possibly arising in this type of flow.
Its double integration using the respective boundary conditions (1d) and (1f) gives the correction of the leading-order velocity profile: