The role of soluble surfactants in the linear stability of two-layer flow in a channel

The linear stability of Couette–Poiseuille flow of two superposed fluid layers in a horizontal channel is considered. The lower fluid layer is populated with surfactants that appear either in the form of monomers or micelles and can also get adsorbed at the interface between the fluids. A mathematical model is formulated which combines the Navier–Stokes equations in each fluid layer, convection–diffusion equations for the concentration of monomers (at the interface and in the bulk fluid) and micelles (in the bulk), together with appropriate coupling conditions at the interface. The primary aim of this study is to investigate when the system is unstable to arbitrary wavelength perturbations, and in particular, to determine the influence of surfactant solubility and/or sorption kinetics on the instability. A linear stability analysis is performed and the growth rates are obtained by solving an eigenvalue problem for Stokes flow, both numerically for disturbances of arbitrary wavelength and analytically using long-wave approximations. It is found that the system is stable when the surfactant is sufficiently soluble in the bulk and if the fluid viscosity ratio $m$ and thickness ratio $n$ satisfy the condition $m<n^{2}$ . On the other hand, the effect of surfactant solubility is found to be destabilising if $m\geqslant n^{2}$ . Both of the aforementioned results are manifested for low bulk concentrations below the critical micelle concentration; however, when the equilibrium bulk concentration is sufficiently high (and above the critical micelle concentration) so that micelles are formed in the bulk fluid, the system is stable if $m<n^{2}$ in all cases examined.

amphiphilic molecule that contains a hydrophilic head and a hydrophobic tail; hence in fluid mixtures such as water-oil systems, surfactants tend to arrange themselves at the interface with their heads in the aqueous phase and their tails in the oil phase. It is also possible for surfactants to be soluble in either of the two phases, either as individual molecules (at sufficiently low concentrations) or colloidal-sized aggregates called micelles (at higher concentrations). The micelles occur only when the concentration is above a critical value called the critical micelle concentration (CMC), and are formed by directing their oil-soluble tails away from the aqueous phase and trapping small droplets of oil inside. Micellisation can therefore promote the mixing of immiscible liquids and can play a crucial role in cleaning processes.
One of the key functions of surfactants is to reduce the surface tension at fluid interfaces. It is known that for bulk concentrations below the CMC, the surface tension decreases with the surfactant concentration according to the Gibbs isotherm (e.g. Chang & Franses 1995). Experimental surface-tension measurements have indicated, however, that the surface tension reaches a saturation level as the surfactant concentration in the bulk becomes larger than the CMC (Song et al. 2006). As soon as the bulk concentration reaches the CMC, any additional surfactant mass is readily available to form micelles and it is expected that the interface and bulk monomer concentrations would remain constant. This implies that at equilibrium and for bulk concentrations above the CMC, there is a disassociation between the change in surfactant mass and the surface tension (which is not sensitive to the concentration of micelles).
The reduction of the surface tension in the presence of surfactant (at sufficiently low concentrations) gives rise to so-called Marangoni forces that act locally on the interface and drive the fluid towards regions of higher surface tension. This phenomenon can significantly influence microscopic multiphase flows where surface tension is the dominant force. It is hence evident that surfactants can be used as a means of manipulating flows and controlling interfacial instabilities (examples of such instabilities arising in film spreading processes are given in the review by ). The ability to control and manipulate multi-layer systems, by either suppressing or enhancing instabilities, is essential in numerous applications that either benefit from a flat and smooth film or require interfacial travelling waves to enhance mass transport; examples include coating applications (Weinstein & Ruschak 2004), oil recovery (Slattery 1974), foam drainage (Shaw 1992) and microfluidics technology ). It is therefore of fundamental importance to understand the physical mechanisms responsible for interfacial instabilities.
The stability of interfacial flows has been studied extensively. In the absence of surfactants, Yih (1967) found that a two-layer parallel flow is unstable to perturbations of long wavelength. The instability occurs as long as the Reynolds number is non-zero and it requires a viscosity contrast across the interface. Yih's work has been extended by other authors to perturbations of arbitrary wavelength (Renardy 1985;Yiantsios & Higgins 1988) and to semi-infinite fluids (Hooper & Boyd 1987) in order to identify the effect of bounding walls. The impact of surfactants on two-layer flows has been analysed predominantly for the case of insoluble surfactants. Linear stability studies indicated that the interface can be destabilised by insoluble surfactants even in the Stokes flow limit (Frenkel & Halpern 2002;Halpern & Frenkel 2003;Blyth & Pozrikidis 2004b;Frenkel, Halpern & Schweiger 2019a,b), and eventually becomes saturated to non-uniform states in the nonlinear regime (Blyth & Pozrikidis 2004b;Pozrikidis 2004a;Wei 2005;Frenkel & Halpern 2006;Bassom, Blyth & Papageorgiou 2010;Samanta 2013;Kalogirou & Papageorgiou 2016;Frenkel & Halpern 2017;Kalogirou 2018).
There have been a number of studies considering the effect of soluble surfactants on the linear stability of falling liquid films (Ji & Setterwall 1994;Karapetsas & Bontozoglou 2013 or two-layer channel flows (Sun & Fahmy 2006;Zaisha et al. 2008;You, Zhang & Zheng 2014;Picardo, Radhakrishna & Pushpavanam 2016). In the latter case, the papers were either based on the simplifying assumption of a non-deformable interface or considered surfactant soluble in both phases, with the surface tension linearly dependent on either bulk concentration. In this paper, we also consider dynamic transport of surfactant at the interface and include the adsorption kinetics. The literature considering flows in the presence of surfactant at high concentrations above the CMC is very limited and is restricted to single fluids with a free surface -see Breward & Howell (2004) for an analysis of steady straining flow of a micellar surfactant solution; Edmonstone, Craster & Matar (2006) for a report on surfactant-induced fingering patterns observed in droplet spreading on clean films; Craster, Matar & Papageorgiou (2009) for a study on break up of surfactant-laden jets; Karapetsas, Matar & Craster (2011) for a review of surfactant-enhanced spreading and superspreading on solid surfaces.
In this article, we present a detailed mathematical model describing a two-layer flow where one of the fluids is contaminated with soluble surfactant. Even though the surfactant could theoretically exist in both fluid layers, here it is only allowed to live in one of the fluids as this scenario is more relevant to practical applications involving water-oil systems, for instance in cleaning processes. The model incorporates surfactant in three different arrangements: molecules adsorbed at the interface, single molecules (monomers) in the bulk fluid and multi-molecule aggregates (micelles) in the bulk. The presence of soluble surfactants introduces additional complexity in the model due to the coupling of the interfacial dynamics with the dynamics in the bulk through mass exchange. To the best of our knowledge, this is the first study that considers a multi-layer flow with soluble surfactants above the CMC. The derived model integrates a number of salient physical properties that play an important role in the dynamics of surfactant-laden multi-layer flows, such as inertia, gravity, surface tension, Marangoni stresses, diffusion and mass transfer between the different surfactant forms. This model reduces to that of Kalogirou (2018) (which is the most general model presented so far and includes the effects of inertia and density stratification) for a two-layer flow with an insoluble surfactant in appropriate limits. We perform a linear stability analysis that yields an eigenvalue problem for the wave speed, and obtain analytical expressions for the two dominant growth rates valid in the long-wave approximation. The linearised system is also solved numerically for perturbations of arbitrary wavelengths and results are presented for situations with bulk concentrations below or above the CMC.
In recent experimental work, Georgantaki, Vlachogiannis & Bontozoglou (2012 investigated liquid film flow with soluble surfactants, and reported the formation of solitary waves preceded by capillary ripples as well as small-amplitude sinusoidal travelling waves. In other related experiments, Shen et al. (2002) studied fibre coating with surfactant solutions and investigated the thickening properties of the film in terms of the bulk surfactant concentration. Hashimoto et al. (2008) considered the flow of surfactant-laden droplets in a microfluidic Hele-Shaw cell and found a range of flow patterns, including elongation of droplets and shear-driven instabilities. However, to our knowledge there are no experiments on two-layer channel flows with surfactants that could serve as a means of validation for our model.
The paper is organised as follows. In § 2 the physical problem is presented and the mathematical model is formulated. Linear stability analysis about the equilibrium state x U Fluid 2: ® 2 , µ 2 Fluid 1: ® 1 , µ 1 H 1 g H H 2 y FIGURE 1. (Colour online) Problem configuration: two superposed fluid layers in a channel of fixed height H, driven by the upper wall motion with velocity U and by a constant pressure gradient G = −(∂p/∂x). The lower fluid is contaminated with surfactant molecules, which can also attach on the interface between the two fluids or form micelles. The size of the molecules is not to scale. is performed in § 3 and asymptotic long-wave expansions are introduced. Numerical results for bulk concentrations below and above the CMC are presented and discussed in § 4. The main conclusions of this study can be found in § 5.

Mathematical model
The flow considered in this paper is illustrated in figure 1. We define a Cartesian coordinate system with horizontal coordinate x and vertical coordinate y, while time will be denoted by t. The problem configuration comprises two immiscible viscous fluids that fill a long horizontal channel, consisting of two impermeable parallel plates at y = 0 and y = H. The two fluid layers are superposed and separated by a distinct interface, each occupying regions 0 y h(x, t) and h(x, t) y H, respectively. The interface at y = h(x, t) can evolve in space and time, but in a steady scenario is flat at y = h 0 H, with 0 < h 0 < 1. The flow is driven in the x-direction by the motion of the upper plate with longitudinal velocity U, as well as a constant pressure gradient G = −(∂p/∂x).

Hydrodynamics
The two fluids have densities ρ 1 , ρ 2 , viscosities µ 1 , µ 2 and thicknesses H 1 , H 2 , where subscripts 1 and 2 refer to the lower and upper fluids, respectively. We define the pressure p j and the velocity vector u j = (u j , v j ) in each region j = 1, 2. The evolution of the flow in each fluid layer is governed by the mass and momentum equations, given by where ∇ = (∂/∂x, ∂/∂y) and g = (0, −g) with gravitational acceleration g.
On both walls the no-slip and no-penetration conditions for the fluid velocities are applied, namely u 1 = (0, 0) at y = 0, and u 2 = (U, 0) at y = H.
(2.2) At the interface the fluid velocities need to be continuous, i.e.
and should satisfy the kinematic condition (2.4) The above kinematic condition is accompanied by a dynamic condition that requires the stress to be continuous across the interface. The stress jump at the interface is given by (Stone & Leal 1990;Milliken, Stone & Leal 1993) where γ is the surface tension and σ j is the stress tensor in fluid j defined by (Batchelor 1967) Also, n is the unit normal vector (pointing into fluid 1), ∇ s is the surface gradient operator and κ is the interfacial curvature, given by (2.5c) Equation (2.5a) takes into account the capillary pressure jump due to interfacial curvature (normal to the interface) and the Marangoni force arising as a result of the variability of surface tension due to the presence of surfactant (in the tangential direction). The stress balance in the normal and tangential directions is provided by taking the dot product of (2.5a) with the unit normal vector n from (2.5c) and unit tangent vector t = (1, h x )/ 1 + h 2 x , respectively, and is equal to where the jump notation f j 2.2. Surfactant transport and connection to surface tension The lower fluid 1 is contaminated with surfactant of bulk concentration C(x, y, t). The surfactant molecules can get transferred onto the interface or form into micelles, whose concentrations are denoted by Γ (x, t) and M(x, y, t), respectively. We suppose that the surfactant molecules join/leave the interface from/to the bulk with adsorption/desorption rates k a , k d , and that the micelles form and break up with kinetic rate constants k b , k m , respectively (figure 2). It is also assumed that the micelles are monodispersed (Shaw 1992) and each micelle has a preferred size consisting of N monomers. We adopt the kinetic model presented by Edmonstone et al. (2006), , which introduces the following fluxes Adsorption k a FIGURE 2. (Colour online) Schematic showing the adsorption and desorption kinetics at the interface, as well as the assembly and disassembly kinetics of micelles in the bulk. The size of the molecules is not to scale.
where the bulk concentration C in J b is evaluated at the interface y = h. The above kinetic model is clearly nonlinear and takes into account that the space on the interface is limited, thus it can become fully packed with monomers at a finite concentration Γ ∞ (the maximum packing concentration), at which no more molecules can be further adsorbed on the interface. The model also assumes that micelles are formed rapidly (since N is typically large) when the critical micelle concentration is exceeded -this is due to the nonlinear term C N . The surfactant monomers at the interface, the monomers in the bulk and the micelles follow appropriate transport equations, given by (Stone & Leal 1990;Karapetsas et al. 2011 with D s , D b and D m denoting the surface, bulk and micelle surfactant diffusivities, respectively. The first term in the interfacial transport equation (2.8) represents the time derivative of a point on the interface which moves in a direction normal to the surface (with an arbitrary tangential componentṙ) -see Wong, Rumschitzki & Maldarelli (1996) for a detailed derivation. Also, u I and u S denote the interfacial velocity and its tangential component, respectively, defined as follows At the bottom channel wall, no-flux conditions are imposed for the bulk and micelle concentrations, i.e. n · ∇C = 0 and n · ∇M = 0, at y = 0, (2.12) and on the fluid 1 side of the interface the following flux conditions are valid (2.13) The first condition in (2.13) describes the exchange of monomers between the interface and the nearby bulk fluid, while the second condition states that the flux of micelles onto the interface in zero. This condition is physically appropriate and essentially means that no micelles can adsorb directly to the interface, but they must first break up into monomers in the bulk. Finally, the total mass M of surfactant in a domain of (arbitrary) length L is conserved, that is (2.14) where micelles can only form as long as the mass available is sufficient for the bulk concentration to exceed the critical micelle concentration.
When the bulk concentration C is below the critical micelle concentration, C CMC , addition of surfactant in the bulk also increases the interface concentration Γ , which in turn reduces the interfacial surface tension γ according to the Gibbs law (Chang & Franses 1995;Pozrikidis 2004b) where R is the ideal gas constant, T is the absolute temperature and C s denotes the bulk concentration at the interface. The Gibbs isotherm assumes an ideal bulk fluid and non-interacting molecules, and is valid when the bulk concentration and the surface concentration are at thermodynamic equilibrium (Hiemenz & Rajagopalan 1997), namely when with γ 0 corresponding to the clean surface tension in the absence of surfactant. Note that even though the above equation of state only relates the surface tension to the interfacial concentration, the surface tension is also affected by the bulk concentration through (2.16). When the concentration in the bulk fluid becomes greater than the critical micelle concentration, then any additional surfactant in the bulk is converted into micelles leaving the interface concentration Γ unchanged and, in turn, the surface-tension constant. The phenomenon of the surface tension settling to a constant value for bulk concentrations beyond the CMC is well known and has been reported in experimental surface-tension measurements -see Elworthy & Mysels (1966), Zhmud, Tiberg & Kizling (2000), Liao, Basaran & Franses (2003), Song et al. (2006). It is important to note that in this argument 'concentration in the bulk' refers to the total concentration C b consisting of both monomer and micelle concentrations, i.e. C b = C + M = C + C N (Danov et al. 1996;Zhmud et al. 2000). Consequently, the plot of the surface tension from (2.17) against C b is seen reaching a plateau at a critical value γ = γ c when C b C CMC (figure 3).

Non-dimensionalisation
The problem is written in non-dimensional form by performing the following transformation of variables where C CMC = (k m /Nk b ) 1/(N−1) (Breward & Howell 2004;Edmonstone et al. 2006;Karapetsas et al. 2011) and L is the length of an (arbitrary) horizontal domain. A number of dimensionless parameters are therefore introduced and are given in table 1. To simplify the notation, we will henceforth drop the superscript stars from all dimensionless variables.

Governing equations within the fluids
The flow in each fluid region is described by the non-dimensional continuity and Navier-Stokes equations, given below for fluids j = 1, 2, The transport of surfactant monomers and micelles in the bulk of fluid 1 is described by the convection-diffusion equations with the dimensionless flux J m given by Fluid property Hydrodynamic Surfactant Solubility ratios parameters parameters parameters The following parameters are also defined: r j = ρ j /ρ 1 and m j = µ j /µ 1 , j = 1, 2, such that r 1 = 1, r 2 = r and m 1 = 1, m 2 = m.
The boundary conditions imposed on the channel walls are The dimensionless forms of the normal and tangential stress jumps at the interface (2.6) are respectively given by The right-hand side terms in both the normal and tangential stress balance (2.26a), (2.26b) indicate the dependence of the surface tension on the local surfactant concentration; the two are connected via the dimensionless nonlinear equation of state with the flux boundary conditions at the interface becoming 29) and the non-dimensional flux J b given by The dimensionless total mass of surfactant is where the domain length L is also scaled in non-dimensional units (by writing L = HL * and then dropping the superscript star from L * ). The model derived in this section reduces to the insoluble surfactant problem by applying one of the limits: (i) Bi → 0, found when there is no desorption from the interface to the bulk, i.e. k d → 0 (Booty & Siegel 2010;Wang, Siegel & Booty 2014). This condition is equivalent to setting the flux J b = 0 in (2.28) and (2.29); (ii) β b 1, which is essentially equivalent to condition (i), see first equation in (2.29) (Jensen & Grotberg 1993;Oron, Davis & Bankoff 1997;Edmonstone et al. 2006;); (iii) K b 1, corresponding to the molecules being attracted to the interface, i.e. k a k d ).
In place of conditions (ii) and (iii), we could instead satisfy R b 1, where R b is a new parameter defined by (2.32) As we will see in later sections, this combined parameter R b will often appear in the stability conditions of this problem. A similar observation has been made by Karapetsas & Bontozoglou (2013) in a related liquid film flow with soluble surfactants.
3. Linear stability analysis 3.1. Equilibrium state At equilibrium both fluid layers have uniform thicknesses, with the interface located at y = h 0 . The basic flow in each fluid layer j = 1, 2 is the standard two-fluid Couette-Poiseuille flow and is given by where p 0 is the undisturbed constant pressure at the interface and the rest of the coefficients are given by Parameter q represents the effect of the pressure gradient. We can define a new parameter for the dimensionless shear rate s by which is equal to the slope of the basic velocity at the interface; we will see later that this shear parameter plays a crucial role in the stability of the problem. The equilibrium state for the surfactant can be found by setting J b = 0, J m = 0 in (2.30) and (2.22), resulting in The total surfactant mass is For any given mass M, we can find the equilibrium concentration Γ by substituting for M and then C from (3.3).

Linearised equations of motion
The interest of this work is to study the linear stability properties of the flow. To formulate the stability problem, we first consider small disturbances (denoted with a hat decoration) to the equilibrium state as follows (3.5a) The incompressible condition (2.19a) allows us to write the fluid velocities in terms of derivatives of scalar streamfunctions, i.e.
The perturbation variables are then written in a normal-mode form where k is the (real) wavenumber, c is the wave speed (generally a complex value) and the tilde variables define the amplitude vector. The amplification or growth rate is therefore given by λ = kIm(c) and hence stability is determined by the sign of Im(c): a negative sign yields exponential decay of perturbations, while a positive sign means linear instability. Linearisation of the kinematic equation (2.25) gives defines the velocity of the moving disturbance relative to the unperturbed interfacial velocity. Equation (3.6) will used in the remaining equations to eliminateh in terms ofΨ 1 (h 0 ). Substituting perturbations (3.5) in the Navier-Stokes equations (2.19b)-(2.19c) and boundary conditions (2.23a)-(2.23b), (after eliminating the pressure) results in the Orr-Sommerfeld boundary value problem for the streamfunction disturbance in each fluid, given by where primes denote differentiation with respect to y. The linearised conditions for vertical and horizontal velocity continuity (2.24) at the interface y = h 0 arẽ Inserting disturbances (3.5) into (2.26) and using the leading-order momentum equations (2.19), yields the linearised normal stress jump at the interface where γ = 1 + β s ln(1 − Γ ) is the equilibrium surface-tension value, and the corresponding tangential stress jump (3.7e) Equation (3.7d) is simplified through cancellation of the second and fourth terms (and by using velocity continuity (3.7c)) when the fluids have equal densities, i.e. for r = 1. Note that with the exception of the right-hand side surfactant term in the tangential stress equation (3.7e), the rest of the linear hydrodynamic system (3.7) up to here is identical to that obtained by Yih (1967) for the stability of multi-layer flow without surfactant.
What is left is to linearise the transport equations for the surfactant. This is achieved by first substituting perturbations (3.5) into equation (2.28), providing the following linear equation for the interfacial surfactant concentration disturbancẽ (3.8) Finally, the convection-diffusion equations (2.20)-(2.21) for the monomers and micelles in the bulk and the relevant boundary conditions (2.23c), (2.29) are linearised; the resulting system for the disturbance in the monomer concentration reads and the corresponding perturbation system for the micelle concentration is (3.10b) Clearly the disturbances for both the monomer and micelle concentrations in the bulk each satisfy a second-order ordinary differential equation with two boundary conditions; while these are homogeneous in the micelle system, one of the conditions (3.9b) in the monomer system is a Robin boundary condition. At this point we should point out that the linear system (3.7)-(3.10) in effect only depends on parameter R b = β b K b and not on β b or K b individually; this can be seen by scalingC → β bC (M needs to be also scaled appropriately) and noting that the factor (K b C + 1) can be written in terms of Γ only (by applying the equilibrium equation (3.3)).

Expansions for long waves
Previous studies on multi-layer flows have shown that the interface is susceptible to long-wave instability, i.e. instability to disturbances with large wavelength. That motivates us to look for a similar instability in multi-layer flows with soluble surfactant, and hence we introduce the following expansions for long waves (i.e. small wavenumber k), c = c 0 + kc 1 + · · · , (3.11a) Ψ j (y) =Ψ j,0 (y) + kΨ j,1 (y) + · · · , (3.11b) Γ =Γ 0 + kΓ 1 + · · · , (3.11c) C(y) =C 0 (y) + kC 1 (y) + · · · , (3.11d) M(y) =M 0 (y) + kM 1 (y) + · · · . (3.11e) The solution procedure followed is to seek solutions to the hydrodynamic and surfactant systems individually at each order until we eventually recover all variables up to order k, and most importantly obtain wave speeds c 0 , c 1 which will provide information of the stability properties of the problem. The leading-order variables in the case of bulk concentrations beyond the CMC are found to be algebraically rather complicated, and therefore the remaining of this section will consider bulk concentrations below the CMC (withM = 0) for the simplicity of presentation.
3.3.1. Hydrodynamic system at O(1) The hydrodynamic system (3.7) at leading order in k is unaware of the presence of surfactant and the solution is exactly the same as the one found by Yih (1967). The leading-order relative velocity of interfacial waves is where the thickness ratio n is connected to the undisturbed lower fluid thickness h 0 through n = (1 − h 0 )/h 0 .

Surfactant system at O(1)
The surfactant equations (3.8)-(3.9) at leading order in k read Combining condition (3.13a) and the second equation in (3.13c) yieldsC 0 (h 0 ) = 0 and (3.13a) can be then used to provide a condition forC(h 0 ) in terms ofΓ 0 . Integrating equation (3.13b) twice in y gives which shows that the leading-order perturbation of the bulk concentration is constant across the fluid. Interestingly, this solution forC 0 is identical to the leading-order bulk disturbance found in liquid film flow down an inclined plate (Karapetsas & Bontozoglou 2014). This is unsurprising considering that the leading-order surfactant system is unaware of the presence of the second fluid.

Surfactant system at O(k)
The bulk disturbance equation and boundary conditions at the next order arẽ Integrating equation (3.15a) twice and using the boundary conditions (3.15b) to determine the constants of integration results in the first-order perturbation for the bulk concentratioñ where the termsΓ 0 ,Γ 1 still need to be found. To obtainΓ 0 we need to consider the interfacial surfactant equation at O(k), and after some algebraic manipulations this is found to be (given for (3.17c) where we have assumed that m = 1. The first factorΓ ins 0 in (3.17a) is identical to the leading-order perturbation in the corresponding insoluble problem, and the second factor S tends to 1 in the insoluble limit R b 1, as expected. The solution forΓ 1 is found by solving the interfacial surfactant equation at O(k 2 ) but it is too lengthy to be given here. Having obtainedΓ 0 andΓ 1 , the leading-order and first-order perturbations for the bulk concentrationC 0 (y) in (3.14) andC 1 (y) in (3.16) are then fully determined.

Hydrodynamic system at O(k)
Analytic calculations for the hydrodynamic system at O(k) are algebraically very cumbersome so the symbolic software Maple is employed in order to find the firstorder perturbation of the wave speed c 1 ; for Re = 0, Bo = 0, q = 0 the solution is the corresponding insoluble component. The solution for c 1 is purely imaginary and hence provides the leading-order expression for the growth rate λ = kIm(c), given by k 2 Im(c 1 ) = c ins Sk 2 . Since the soluble component S tends to 1 in the limit R b 1, the growth rate reduces to the insoluble rate. As mentioned before, the expression (3.18a) for the growth rate only depends on parameter R b = β b K b through S given in (3.17c) and not on parameters β b or K b individually. It is important to note that expressions (3.17) and (3.18) are only valid for m = 1 and (provided that m > 1). At m = 1 it was shown by Frenkel & Halpern (2002) and Halpern & Frenkel (2003) for the corresponding insoluble problem that the long-wave expansion for the wave speed is not regular, but the first-order term in the expansion is of O(k 1/2 ) instead of O(k). This happens because in this particular limit, the linear term in the quadratic equation that determines the wave speed c does not contribute to the leading-order solution in k. We have proved that a similar result holds when R b = R * b , in which case the first-order perturbation of the wave speed is given by 3.3.5. Second (Marangoni) mode The mode found above is a 'hydrodynamic mode', in the sense that it exists independently of the presence of surfactant; this mode is known to be stable in the Stokes flow limit and in the absence of surfactant (Yih 1967) but it is generally unstable for non-zero Reynolds number. Several authors in the past have reported a second 'Marangoni mode', associated with the surfactant concentration at the interface and coming from the interfacial transport equation (e.g. Frenkel & Halpern 2002;Pereira & Kalliadasis 2008;Karapetsas & Bontozoglou 2014). We can find the corresponding Marangoni mode for this problem by expanding the wave speed as in (3.11a) and introducing expansions for the streamfunctionsΨ j with a leading term of order α −1 (we note that the expansion is the same as in the insoluble problem, see Frenkel & Halpern (2002)), while the surfactant concentrationsΓ andC are expanded with a leading term of order α −2 . We obtain the following solution for the wave speed at leading order whereas the expression of the O(k) term c 1 is rather unwieldy and will be omitted. Clearly, the leading-order wave speed c 0 is real and therefore the growth rate λ = kIm(c) passes through k = 0. In the insoluble limit R b 1, the above expression c 0 → −qh 2 0 + sh 0 = qh 2 0 + wh 0 (upon setting n = (1 − h 0 )/h 0 and using (3.2)) which is equal to the undisturbed interfacial velocity u 1 (h 0 ).
The second mode has a similar form as (3.18); it also includes the solubility factor S (defined in (3.17c)) and exhibits a different asymptotic behaviour at critical point R * b , similar to the first mode presented in § 3.3. Here we discuss the behaviour of these two modes as the viscosity ratio m varies. In figure 4, the two long-wave modes are shown with solid lines and the corresponding insoluble modes are shown with dashdotted lines -note that the vertical asymptotes do not signify a breakup of the longwave theory but rather a change in the asymptotic expansion as discussed above. Both modes are seen to change stability as the boundary m = m * > 1 is crossed, but the overall system always remains unstable in the vicinity of m * (similar to the behaviour of the corresponding insoluble system at m = 1). The critical point m * exists as long as R b > 1/(1 − Γ ) 2 n 2 and tends to 1 in the insoluble limit R b 1.

Other normal modes
The analogous insoluble problem admits only two normal modes under conditions of Stokes flow. Both of these modes vanish when the disturbance has infinite  wavelength, that is when k = 0. When the surfactant is soluble there is an infinite number of modes that originate from the bulk transport equation (3.9) -see appendix A for a method of obtaining these modes in the case of Stokes flow. Two modes exactly pass through the origin k = 0 but all the other modes reach a constant (negative) value at k = 0; further analysis of these remaining modes is provided in appendix B.

Numerical method
A general solution of the linear eigenvalue problem (3.7)-(3.10) for perturbations of arbitrary wavelengths can only be obtained numerically. Here we use the Chebyshev collocation method which requires to first map each fluid region to the canonical interval [−1, 1]. This is achieved by performing a y-coordinate transformation to new coordinates y 1 , y 2 ∈ [−1, 1] defined by with the interface located at y 1 = y 2 = −1. The streamfunctionsΨ 1 ,Ψ 2 and surfactant concentrationsC,M are then written in terms of Chebyshev expansions (Orzag 1971;Boomkamp et al. 1997) where T (y j ) is the th Chebyshev polynomial of the first kind and N 1 , N 2 are appropriately chosen truncation levels for each fluid layer j = 1, 2. The total number of unknowns considered is 3N 1 + N 2 + 5, comprising N 1 + 1 expansion coefficients ψ 1 , N 2 + 1 coefficientsψ 2 , N 1 + 1 coefficientsc , N 1 + 1 coefficientsm andΓ . The linearised set of equations (3.7)-(3.10) is assembled into a system (A · x) = c(B · x), where A, B are square matrices of size 3N 1 + N 2 + 5. The linear system then becomes a finite-dimensional generalised eigenvalue problem for the complex eigenvalue c = c r + ic i and the associated eigenvector. The solution of this system provides the growth rate λ = kc i for any wavenumber k.
In practice, the two Orr-Sommerfeld equations are solved with N 1 − 3 and N 2 − 3 interior collocation points, respectively, while the bulk and micelle concentration equations are solved at N 1 − 1 interior points each. The remaining 13 equations come from the boundary conditions (3.7b)-(3.7e), (3.8), (3.9b), (3.10b). Note that some of the boundary conditions do not contain the eigenvalue c and so a number of rows in B are zero, making the matrix singular. The singularity in B in general leads to the emergence of spurious eigenvalues (Dongarra, Straughan & Walker 1996), which were carefully eliminated in order to obtain accurate numerical results. This was achieved in one of the following ways: (i) by removing the null rows and columns and obtaining a reduced system; or (ii) by using the built-in generalised eigenvalue solver eig in MATLAB, which implements a QZ algorithm. Computations with both of these methods have been carried out with eigenvalues being in excellent agreement in both cases.
The system at hand has an infinite number of eigenvalues and associated eigenfunctions for non-zero values of the Reynolds number. In the Stokes flow approximation, Re = 0, the corresponding problem with insoluble surfactant has precisely two eigenvalues. As already mentioned, the problem with soluble surfactant considered here has infinitely many eigenvalues even for Stokes flow; the two dominant modes are identified and demonstrated in the numerical results that follow in this section.
To validate the code and check its numerical accuracy, a number of test cases were simulated and numerical solutions were compared to available results from the literature. As a first test case, growth rates were obtained for multi-layer Couette or Poiseuille flow with an interface that is devoid of surfactants, and comparisons with results presented in Renardy (1985) and Yiantsios & Higgins (1988) showed excellent agreement. A further check was performed by computing growth rates for single-layer or multi-layer flow with insoluble surfactant at the interface and recovering those obtained by Halpern & Frenkel (2003), Blyth & Pozrikidis (2004a) and Pereira & Kalliadasis (2008). We have also solved a modified code for liquid film flow in the presence of soluble surfactant and growth rates were confirmed to be identical to those found by Karapetsas & Bontozoglou (2013). Furthermore, in the case of multi-layer flow with soluble surfactant it was verified that the code reproduces the results of an independently written code for the corresponding insoluble problem when Bi = 0, or Bi 1 and R b 1.
As a final validation case we compared the numerically computed growth rates to the long-wave asymptotic solutions presented in § 3.3. Figure 5 demonstrates an example, where the two dominant modes calculated using the Chebyshev method are shown with solid black lines and the long-wave predictions are portrayed with red circles. Clearly the numerical and asymptotic growth rate curves become indistinguishable as they approach the origin. We verified that this is indeed the case for all other results presented next.

Numerical results
The primary aim of this section is to identify the effect of surfactant solubility and sorption kinetics on the stability of the interface, in order to distinguish the stability properties of multi-layer flow with soluble surfactant from those of the equivalent flow where the surfactant is treated as insoluble. The degree of surfactant solubility is represented in the mathematical model by parameter R b , with R b 1 being highly soluble in the bulk and R b 1 nearly insoluble (we have confirmed that β b or K b do not affect the numerical results individually but only through the parameter R b = β b K b , as expected). A measure of sorption kinetics is provided by the Biot number Bi, with Bi 1 corresponding to surfactant leaving the interface very slowly and Bi 1 indicating a fast desorption rate towards the bulk fluid. We will investigate the influence of the two parameters Bi, R b on the stability properties of the system and explore how these properties are affected when the bulk concentration becomes larger than the CMC.
Considering the large number of dimensionless parameters involved in the problem at hand (table 1), it is practical to fix some of these parameters in the numerical simulations. All the results presented in this section will be for Stokes flow, Re = 0, and unless otherwise specified the following parameters will be kept fixed: r = 1, Bo = 0, Ca = 0.1, Ma = 0.1, Pe s = 1 × 10 8 , Pe b = 100, Pe m = 100, K m = 1 (note that the last two parameters are relevant only in the presence of micelles). Frenkel & Halpern (2002) showed that the interface between two viscous fluids under Stokes flow conditions exhibits long-wave instability in the presence of insoluble surfactant, but only if an underlying shear is imposed. We investigated if the related problem with soluble surfactant follows the same behaviour and indeed we found that the problem is stable to small disturbances of any wavelength if the basic flow shear is removed. This result was shown using both numerical calculations  of the linearised system and long-wave analysis with shear rate s = 0 (results not shown for brevity). It is therefore imperative that an underlying shear flow is present for this problem to exhibit any kind of instability to interfacial disturbances. The rest of the results presented in this section will accordingly employ a non-zero shear, s = 0, given by (3.2). We start the numerical investigation by considering a parameter set that supports instability when the surfactant is mostly attracted to the interface (i.e. insoluble or nearly insoluble) and slowly strengthen the effect of surfactant solubility or sorption kinetics by varying parameters R b or Bi. Figures 6 and 7(a) use m = 0.5, n = 10 (h 0 = 1/11) and Γ = 0.5, which corresponds to an unstable interface in the case of insoluble surfactant (since m < n 2 , see Frenkel & Halpern (2002)). In figure 6, the effect of solubility parameter R b on the stability of the interface is displayed. For small and fixed Bi = 0.01 (figure 6a), the dominant growth rate is close to the one for the insoluble problem when R b = 100 (shown with a thick solid line); decreasing R b = 10, 1, 0.2, 0.1 reduces the cutoff wavenumber and eventually stabilises the problem. The growth rates are seen to pass through the origin and to have a non-monotonic behaviour even after they are stabilised, with instability established for sufficiently long waves only. For larger Bi = 1, the stabilisation occurs at a higher value of R b and is monotonic as R b decreases (figure 6b shows the two most significant modes for each value of R b = 10, 2, 1, 0.5). When R b becomes smaller than the value R b = 0.25 (not shown in the figure), the two curves cross each other and the second mode becomes dominant (but both modes remain stable). Fixing R b = 5 and increasing Bi = 0.001, 0.0025, 0.005, 0.01 reduces the range of unstable wavenumbers (but only by a small amount) as shown in figure 7(a), whereas larger values of the Biot number leave the growth rate curves unaffected. In the large Biot number limit the transport is controlled by diffusion and the interface-bulk exchange kinetics are in equilibrium (Booty & Siegel 2010).

Bulk concentrations below the CMC
A particularly interesting result is seen when varying the equilibrium surfactant concentration at the interface, Γ . In the insoluble problem, increasing the basic concentration Γ would make the problem more unstable by increasing the cutoff wavenumber. Here this is not the case, as illustrated in figure 7(b) for fixed R b = 1, Bi = 1 and Γ = 0.1, 0.3, 0.5, 0.7, 0.9 (again, m = 0.5, n = 10); clearly the growth rates follow a non-monotonic behaviour as Γ is varied, with instability supported only for Γ < 0.723. An increase in the basic concentration Γ enhances the range of instability as long as the concentration is small and between 0 < Γ < 0.42. Further increase in Γ beyond the value 0.42 is seen to diminish the interval of unstable wavenumbers and leads to complete stabilisation at Γ = 0.723. A similar behaviour has been reported in the study of falling film flow in the presence soluble surfactant, see Karapetsas & Bontozoglou (2013). The next set of results demonstrate that surfactant solubility can destabilise a system which is stable for insoluble surfactant. The flow with insoluble surfactant is known to be stable when m n 2 (Frenkel & Halpern 2002); computations are hence performed for m = 3 and n = 1.5 (i.e. h 0 = 0.4), so as to satisfy this condition. Parameter Γ is fixed to the value 0.5 while C is found by (3.3). Figure 8 illustrates the growth rates for fixed R b = 1 and increasing Bi = 0.001, 0.01, 0.025, 0.1, 1. The thick solid line for Bi = 0.001 is nearly identical to the insoluble growth rate and is stable, but as the Biot number increases the growth rates become unstable. For larger values of the Biot number the growth rate curves are found to be nearly indistinguishable from the curve for Bi = 1. Results for varied R b are considered next and we fix Bi = 1. Figure 9(a) depicts the dominant growth rate for R b = 100, 10, 3, 2, 1. Lowering R b from a large value R b = 100 decreases the stable dominant mode and at the same time increases the second most dangerous mode which is also stable (not shown in the figure). This behaviour persists until around R c b = 2.4426, at which point the two growth rates become almost identical in a small region near the origin but both are stable (although the dominant growth rate curve has an inflection point). The dominant mode eventually crosses the k-axis and becomes unstable at the critical value R s b = 2.4312. These results are in line with the predictions of the long-wave theory, as can be seen in figure 9(b): there are two negative modes for R b > R s b , which are seen to cross each other at R c b = 2.4426 and the second mode is positive between R  where R * b = 2.42 (as given by (3.19)) is shown with a vertical asymptote in the figure. At the value R * b = 2.42 the long-wave expansion has a different form, according to the discussion in § 3.3.4. The system is unstable for all 0 < R b < R * b , but as R b → 0 the positive growth rate tends to 0 from above, i.e. it becomes almost neutral.
All the results presented so far considered some parametric sets which support instability for sufficiently long waves. However it has been found by Halpern & Frenkel (2003) that when the surfactant is insoluble, it is possible for instability to exist in a finite interval of wavenumbers bounded below away from the origin  (the so-called mid-wave instability), while long and short waves are stable. We find a similar result here for soluble surfactant (note that the mid-wave instability has also been reported for related systems in Picardo et al. (2016) and Frenkel et al. (2019b)). In figure 10 we take m = 17, n = 4 (i.e. h 0 = 0.2) and Γ = 0.5; for these parameter values, the insoluble problem is long wave stable since m n 2 but exhibits mid-wave instability. The numerically computed growth rates for fixed R b = 0.1 and increasing Bi = 0.001, 0.01, 0.1, 1, 10 are shown in figure 10(a). The range of unstable wavenumbers reduces as the Biot number increases but only until Bi = 0.526, above which the instability interval vanishes and the system is stabilised. Keeping the Biot number fixed to Bi = 0.1 while lowering R b = 100, 10, 1, 0.1, 0.01 reduces the length of the unstable wavenumber interval, until around R b = 0.081 when long-wave instability also arises -see figure 10(b) and the inset therein. A similar study but now for larger Bi = 1 with decreasing R b = 100, 1, 0.5, 0.1, 0.05 depicts that the range of wavenumbers which support instability is reduced until the critical value R b = 0.126, below which the interface becomes completely stable (figure 10c). The problem remains stable as far as R b = 0.081, at which value instability for long waves emerges (note that this critical value is the same for both values of Bi shown here as the mode that becomes unstable is independent of the Biot number (cf. § 3.3.4)).

Bulk concentrations above the CMC
If the surfactant concentration in the bulk exceeds the CMC, micelles are formed in the lower fluid. In the case where the interface is not subject to any shear, the system is found to be stable to perturbations of any wavelength. For non-zero shear rate, s = 0, we conducted numerical experiments aiming to determine the effect of micelle formation on the stability of the flow. This is achieved by changing the amount of available surfactant mass at the equilibrium state, M, starting with bulk concentrations below the CMC (for small M) and gradually increasing the concentration until C > 1 so that the bulk monomer concentration is above the CMC. Given a value for the total mass M, we find from (3.3), (3.4) that Γ satisfies a polynomial of degree N + 1, with only one physically acceptable real root in the range 0 < Γ < 1, so that the equilibrium state is unique. Figure 11 shows the dominant growth rates for various values of the surfactant mass, M, and for two different micelle sizes N = 10 and N = 50 (we also set m = 0.5, n = 10). While the available mass is relatively small, we find that the growth rate curves are indistinguishable from the corresponding curves found in the absence of micelles; this is physically anticipated as at low mass availability the bulk concentration is small and below the CMC. For the chosen parametric set and N = 10, this is true for approximately M 0.3. Increasing the mass diminishes the range of unstable wavenumbers and stabilises the system for M 0.519 (figure 11a). When the preferred micelle size is larger at N = 50, more mass needs to be available to stabilise the interface, as shown in figure 11(b). The growth rate remains identical to the equivalent soluble one (obtained forM = 0) up to M = 0.4, but a rather rapid transition to stabilisation is seen as the surfactant mass increases to M = 0.57. The final stable growth rate for large M is confirmed to be identical to the growth rate for a clean system without surfactant.
The latter observation is confirmed in figure 12 which illustrates the variation of surfactant concentrations Γ , C, M at equilibrium as functions of the total mass M, for micelle size N = 50. It is noteworthy that micelles only start to form for mass higher than approximately M = 0.57, after which both C and Γ saturate to constant values while M keeps increasing (the saturated values are approximately C → 1 and Γ → 0.5). This implies that as soon as the micellisation starts, the adsorption and desorption processes are suspended, leaving the interface at a constant surface concentration while any additional mass causes the micelle concentration to grow (Burlatsky et al. 2013). The stabilising action of micelles is physically associated with the plateau in the interfacial tension for bulk concentrations beyond the CMC (see figure 3), consequently the multi-layer Stokes flow is expected to be stable (similarly to the surfactant-free problem with constant surface tension). In cases with m n 2 , it is observed that the formation of micelles reduces the range of unstable wavenumbers and drives the growth rates towards zero; this is found to be true even for large values of mass M or micelle size N (data not shown), and was corroborated using the long-wave asymptotic results.

Conclusions
We have presented a mathematical model for the dynamics of a two-layer surfactant-laden flow that allows for surfactant present at high concentrations (above or below the CMC). The model comprises governing equations for the hydrodynamics and appropriate transport equations for the surfactant concentration at the interface, the concentration of monomers in the bulk fluid and the micelle concentration. It accounts for all the salient physical effects, including inertia, density stratification, viscosity contrast, Marangoni stresses, surface and bulk diffusion, adsorption and desorption kinetics and micellar dis/assembly kinetics. In the limit of vanishing desorption or rapid adsorption rates, the model reduces to that of Kalogirou (2018) for a two-layer flow with an insoluble surfactant.
We have performed a linear stability analysis and solved the resulting Orr-Sommerfeld eigenvalue problem to determine the growth rate. This was done analytically in the limit of long-wave disturbances (assuming surfactant concentrations below the CMC) and numerically for perturbations of arbitrary wavelengths using a Chebyshev collocation method. We note that while the Orr-Sommerfeld system includes inertia and gravitational effects, in our results we excluded both of these. The numerical investigation focused on the effect of surfactant solubility and/or sorption kinetics on the stability of the interface by varying the two key parameters: R b , which is the ratio of the adsorption rate to the desorption rate; or Bi, which is proportional to the rate of desorption.
We have presented numerical results for Stokes flow and for bulk concentrations below or above the CMC. In both cases, short waves were seen to be stable as expected due to the action of surface tension forces (Hooper & Boyd 1987), but disturbances of long or intermediate wavelengths were found to be unstable under certain conditions; we have referred to these as long-wave and mid-wave instabilities, respectively. The instability due to a soluble surfactant is manifested provided that a non-zero shear rate is maintained at the interface. We have analysed the effect of sorption kinetics or surfactant solubility by either increasing the Biot number Bi or reducing the value of R b , respectively. It was found that a flow that is long wave unstable in the presence of insoluble surfactant may be stable if the surfactant is soluble (of concentration below the CMC). Karapetsas & Bontozoglou (2014) suggested that the stabilisation due to surfactant solubility is connected to a near 90 • phase shift between the interface displacement and the mass flux J b (see their figure 3); this has the effect of redistributing surfactant at the interface so as to stabilise the system. We have confirmed that a similar phase shift is observed here. In particular, the maxima in the eigenfunction of the interfacial disturbance occur at almost the same spatial locations as the maxima of the interfacial concentration disturbance, and this mitigates the Marangoni forces. In the case where the corresponding flow with insoluble surfactant is stable (m n 2 ), the role of soluble surfactant was found to be destabilising. In this case, the maxima in the interface displacement and the minima in the interfacial concentration disturbance almost coincide, leading to the intensification of the local Marangoni stresses. Interestingly, long-wave instability was found to co-exist with mid-wave instability at specific parametric sets.
When the lower fluid is not the most viscous fluid, i.e. for m 1, the long-wave approximation of the dominant growth rate has been seen to change asymptotic behaviour from k 2 to k 3/2 at certain values of the parameters (such as R b = R * b or m = m * ). This was confirmed both numerically and analytically and is similar to what Halpern & Frenkel (2003) found for multi-layer flows with insoluble surfactant.
The solubility factor S in (3.17c) is less than one for m < 1, causing the reduction of the leading-order interfacial concentrationΓ 0 and hence weakening the action of the Marangoni forces. Factor S becomes smaller (but always remains positive) as R b is reduced, therefore the attenuation of the Marangoni stresses has a stabilising influence as the solubility effects become stronger. We also note that as the critical points R b = R * b or m = m * are crossed, factor S changes sign, resulting in a 180 • change in phase of the interface concentrationΓ 0 ; the two dominant modes also change sign (and phase), thus the system is always unstable in the vicinity of the boundary points (cf. figures 4 and 9b). For R b > R * b or m > m * , the solubility factor S is greater than one, and therefore the Marangoni effects are strengthened with increasing solubility.
When the total mass of surfactant is small, the surfactant concentration in the bulk remains below the CMC. For higher values of available mass, the bulk concentration eventually becomes larger than the CMC, micelles start to form and their action is seen to stabilise the system very rapidly. Any additional mass only increases the concentration of micelles while the interface and bulk concentrations saturate to constant values. Consequently at surfactant concentrations beyond the CMC, the interface is virtually at a constant surface tension, the Marangoni stresses are eliminated and the system behaves as if it was clean (and in the case of Stokes flow, it is stable).
This work represents the first attempt to examine the stability of a multi-layer flow with soluble surfactants beyond the CMC, by incorporating the presence of micelles as well as their formation and break up. The proposed mathematical model provides a useful framework for modelling the nonlinear dynamics of interfaces in multi-layer flows with surfactant above the CMC. This study can motivate experimental work to investigate the predictions of the linear stability theory, as well as nonlinear simulations to examine the dynamics beyond the linear regime. Such nonlinear investigations are currently in progress.