Transport of soluble surfactant on and within a foam ﬁlm in the context of a foam fractionation process

This study models soluble surfactant transport on and within a foam ﬁlm during a foam fractionation process. Marangoni ﬂow drives surfactant onto the ﬁlm, and also surfactant concentration is assumed to be uniform across the thin foam ﬁlm. Adsorption isotherms are used to couple mass transfer equations, so as to determine the evolution of the total amount of surfactant (surface plus bulk) at any ﬁlm location. Surfactant transport is considered both in the absence and presence of ﬁlm drainage. It is observed that having soluble surfactant slows down evolution of Marangoni-driven ﬂow compared to cases assuming insoluble surfactant. This is because in soluble surfactant cases, surfactants diffuse to the bulk even whilst being transferred by Marangoni ﬂow onto the ﬁlm surface. Furthermore, it is observed that a quasisteady state typically occurs after a long time, such that Marangoni ﬂow and ﬁlm drainage ﬂow become comparable.


Introduction
The separation of surface-active components is essential in numerous industries, such as waste treatment, food, pharmaceutical, environmental-related, and many more (Chang et al., 1992;Gerken et al., 2006;Matsuoka et al., 2022). One of the important methods to separate surfactants from an aqueous solution is foam fractionation (Shedlovsky, 1948;Lemlich, 1972;Keshavarzi et al., 2022). This physicochemical process is often competitive compared to other methods used in this area, namely ultrafiltration, ion exchange, precipitation, membrane filtration and coagulation (Cheang and Zydney, 2003;Schmitt, 2001). This is primarily due to the simplicity of the equipment, low cost and mild operating conditions, as well as environmental compatibility (Du et al., 2000;Mukhopadhyay et al., 2010). Furthermore, the applicability of this method in the separation of dilute solutions, beyond the limits of the other techniques, has made it an attractive alternative (Grieves, 1975;Rubin and Melech, 1972). Indeed the foam fractionation process has various applications, such as radioactive effluent purification (Bergeron and Walstra, 2005), separation of hydrophobic proteins and enzymes (Keshavarzi et al., 2022;Linke et al., 2007;Ackermann et al., 2003;Crofcheck et al., 2003), separation of non-polar compounds (Backleh-Sohrt et al., 2005), production of pharmaceutical products (Lockwood et al., 1997;Datta et al., 2015), environmental problem remediation, (Buckley et al., 2022;Jones and Robinson, 1974) and food processing (Chang and Franses, 1995).
Foam fractionation is based on the adsorption of a surfactant on the bubble surface, as bubbles rise through a solution (Lemlich, 1968). The recovered liquid from the top of the column (foamate) is richer in surfactant than the remaining liquid (Stevenson and Jameson, 2007). Moreover, utilization of reflux can increase significantly the separation ability of a foam fractionation device, and thus increase the concentration of the foamate (Lemlich and Lavi, 1961;Martin et al., 2010;Stevenson et al., 2008;de Lucena et al., 1996). In a fractionation column with reflux, part of the surfactant rich foamate is collapsed and returned to the column. Then, it drains through Plateau borders, where it mixes with any liquid already in those Plateau borders and increases the Plateau border surfactant concentration. This in turn can enhance the surfactant concentration on the foam film surfaces and within the bulk of the foam films. As a result, reflux improves the overall foam fractionation process efficiency. This then is the reason why in this study, we focus on the foam fractionation process with reflux.
This study considers the case of comparatively dry foams with typically polyhedral bubbles (Matzke, 1946;Exerowa and Kruglyakov, 1997). Foam films can then be treated as being thin (i.e. thickness much smaller than length) and are in contact at their edges with Plateau borders. Fig. 1 shows a schematic figure of such a film, indicating how either insoluble or soluble surfactant might be distributed on and within the foam film. Here C and C Pb are surfactant surface concentrations on the film and in the Plateau border, and c and c Pb are analogous bulk concentrations. Due to reflux, C is typically different from C Pb , and c is typically different from c Pb .
There are various physical mechanisms influencing surfactant transport from Plateau border to film and vice versa. For instance, Marangoni flows due to differences between C and C Pb , leading in turn to differences in surface tension, have been studied by numerous authors (Vitasari et al., 2013;Karakashev et al., 2010;Elfring et al., 2016;Vitasari et al., 2016;Grassia et al., 2016). Several studies have also been carried out on film drainage and its effects on foam and foam film stabilization (Coons et al., 2003;Karakashev and Nguyen, 2007;Breward and Howell, 2002). By the same token, the transport of surfactant onto the foam film in the presence of both Marangoni forces and film drainage has been studied (Vitasari et al., 2013;Karakashev et al., 2010;Elfring et al., 2016;Vitasari et al., 2016). However, even though we know surfactants are soluble (Li et al., 1998;Stevenson, 2012), in most previous studies, surfactant transport into the bulk of the film has been neglected (Jensen and Grotberg, 1992;Bruell and Granero-Belinchón, 2020;Vitasari et al., 2013;Jensen, 1995;Yeo et al., 2003). Many typical surfactants have nevertheless a substantial solubility in water which can then impact on flow behaviour (Roché et al., 2014;Tanaka et al., 1973;Schmitt, 2001;Fainerman et al., 2001).
In the case of insoluble surfactants, the work of Vitasari et al. (2013) developed a model for how quickly surfactants move onto the foam film, which is schematically presented in Fig. 1a. This model has been developed both in the absence and the presence of film drainage and has revealed that the surfactants can usually move onto the foam films quickly compared to typical bubble residence times in a fractionation column. However, when film drainage is present, it competes against the Marangoni effect and slows down the surfactant transport on the foam surface. An equilibration of surfactant concentration between the Plateau border and the film is achieved after a sufficiently long time, at least in the absence of film drainage. Moreover, in the presence of film drainage, a quasisteady-state surfactant concentration is reached with a lower surfactant concentration in the film than in the Plateau border. Despite all the achievements of this model, the lack of consideration of any surfactant inside the bulk of the film has led us to present the current model, in which surfactant solubility and its transport inside the bulk are considered (Fig. 1b). In addition, a model has been developed here for the behaviour of the soluble surfactant approaching a quasisteady-state configuration at long times, albeit compared to a model already presented by Vitasari et al. (2013), some modifications are required in order to account for surfactant solubility.
The present model considers in particular the diffusive transport of soluble surfactant into the bulk of the foam film. Diffusive transport from surface to bulk can in principle occur at different rates relative to Marangoni and film drainage flow, but the specific limit that we will look at here is when surfactants can diffuse very quickly across a thin film, even though any diffusion along the film is very slow. As we explain later on, mathematically, we can express this in terms of the product of a Péclet (Pe) number and a film aspect ratio (D), with this product being a small parameter. Hence, we call it a small Pe D limit. To use the model to investigate Here C and C Pb are surface concentrations of surfactant on the film and Plateau border, while c and c Pb are the bulk concentration in the film and Plateau border, respectively. In the present work, C on the film depends on coordinate x and time t. Also, c in the film depends on x and t, but is taken to be almost independent of coordinate z. The figure is not to scale as in reality, the film thickness is much smaller than its length. the performance of foam fractionation for a soluble surfactant, we define (later on) recovery and enrichment factors which are quantities often measured in fractionation studies (Gerken et al., 2006;de Lucena et al., 1996;Stevenson, 2014;Stevenson, 2012).
This study is laid out as follows. Firstly, Section 2 deals with the theory used in mathematical models of soluble surfactant transport. That section is divided into several subsections and explains important dimensionless parameters and governing equations and uses them in the so-called small Pe D limit to develop an equation for the evolution of the total amount of surfactant at every position along a film. After that Section 3 deals with the numerical simulation approach to solve the model and explains some benchmarking methods to check the validity of the simulation approach. Then, Section 4 presents results and discussion, while Section 5 provides the conclusions of this study.

Mathematical model for soluble surfactant transport
In this section, we revisit the relevant transport equations that have been developed by Vitasari et al. (2013), albeit extending that work to consider the possibility of surfactant being present in the bulk. Note however that the work of Vitasari et al. (2013), whilst retaining much of the essential physics of flow and mass transfer on foam films, involves a number of significant simplifications, and we will highlight these as we develop the model to be used here.
Section 2.1 identifies key dimensionless groups controlling the behaviour of the system under consideration. Here we utilise dimensionless groups previously identified by Vitasari et al. (2013). However, a newly introduced solubility parameter is a key parameter also employed in our model. Then, we look at the flow fields in Section 2.2. Since the physical mechanisms that drive the flow fields are Marangoni flow and film drainage, some background to these mechanisms is discussed in the supplementary material Secs. S 1 and S 2. Section 2.3 describes mass transfer equations on the foam film and within the bulk. Modelling of insoluble surfactant transport on the foam surface only is reviewed in supplementary Sec. S 3. However, returning to the soluble case in Section 2.4, analysis of the small Pe D limit, which is the main novel contribution of this study compared to previous studies, is presented. Section 2.5, meanwhile, consists of a new approach to obtaining a linear adsorption isotherm by modifying the Langmuir adsorption isotherm. In particular, we start with a Langmuir isotherm and make approximations, which we call respectively a global Henry isotherm and local Henry isotherm. This approach enables us to have the convenience of working with a linear isotherm, albeit one which is still a good representation of the original Langmuir isotherm over a concentration range of interest. Detailed explanations of these approximations can be found in supplementary Sec. S 4. Subsequently, Section 2.6 presents a new variable called ''total amount of surfactant at any film location", which accounts for surface and bulk taken together at a specific position along the foam film. This then is the variable for which we solve. The solutions are ultimately obtained in terms of variables recast in dimensionless form, as described in supplementary Sec. S 5 with a solution being obtained numerically via a method outlined in supplementary Sec. S 6. Parameters to use within the solution are discussed in Secs. S 7 and S 8. Another key quantity that accounts for the overall amount of surfactant on the surface and within the bulk in a specific time has been introduced in Section 2.7: this quantity is measured over the entire film, not just at a specific location. Moreover, surfactant effective concentration is presented in Section 2.8, giving a measure of the extent to which adsorption enhances surfactant concentration over and above bulk liquid. This then leads within Section 2.9 into a discussion of sur-factant recovery and enrichment during fractionation. Finally, Section 2.10 deals with the system's behaviour at late times, while the details of those late-time solutions can be found in supplementary Sec. S 9.

Dimensionless groups
Here key dimensionless groups that appear in the model are presented. One of the parameters controlling the behaviour of the system is the Péclet number (Bird et al., 2007) which is obtained here based on comparing convective surfactant transport by the Marangoni effect to the diffusive transport of surfactant (see supplementary Sec. S 1.1). In the present model, the Péclet number (Pe) is expressed as follows: where G is Gibbs elasticity parameter (surface tension change per change in logarithm of surfactant surface concentration), d 0 is initial film half-thickness, l is liquid viscosity, L is film half-length and D is diffusion coefficient. Péclet number can be thought of as being the ratio between characteristic velocity for Marangoni flow (along the film) and characteristic velocity for diffusion (across the film) and is typically a relatively large number (see Table S 2). Meanwhile, the initial aspect ratio between film half-thickness and film halflength is defined as: where d 0 is initial film half-thickness, and L is film half-length. Here D is typically a very small number as the film is considered to be extremely thin. This small number then causes the product of Pe and D to be a relatively small number. The other significant parameter is the solubility parameter which is defined as below: where d 0 is initial film half-thickness, C Pb is surfactant surface concentration at the Plateau border, and c Pb is surfactant bulk concentration at the Plateau border. This is a key parameter for this work, describing the amount of dissolved surfactant relative to the amount of surfactant on the surface. This parameter is physically the film thickness relative to a depletion length associated with transferring surfactant between bulk and surface. Insoluble surfactant has S ! 0, whereas a highly soluble one has S ! 1. Note that for values of S ( 1, even though it is possible to compute the amount of surfactant in the interior of the film, there is so little surfactant in the interior compared to the surface that one might as well treat the surfactant as totally insoluble. By contrast, computing what is happening in the interior of the film is much more relevant for larger S values. For the parameters that we examine here, it turns out that S is an order of magnitude or so greater than unity (see Table S 2).
As already alluded to, film drainage can also be relevant to the surfactant transport process (see supplementary Sec. S 2). Hence, we also define a dimensionless film drainage velocity parameter, denoted V R . It is defined such that the typical ratio between drainage flow velocity (Dowson, 1987) and Marangoni flow velocity is on the order of V R . Here V R is expressed as follows: where d 0 is the initial film's half-thickness, c Pb is surface tension at the Plateau border, a is the Plateau border's radius of curvature and G is Gibbs elasticity. This is usually a small parameter (see Table S 2), suggesting that film drainage-driven surfactant transport is typ-ically rather slower than Marangoni-driven surfactant transport. That said, in the case of soluble surfactant, it turns out that there are mechanisms by which Marangoni-driven transport can be slowed down also, so that in relative terms, film drainage-driven transport might become important even though V R is small. This is a point to which we return later on. Remember also that the focus here is upon foam fractionation with reflux. Increasing the amount of reflux flowing through a fractionation column causes the Plateau borders to swell, increasing a and hence reducing V R and thus reducing the impact that any film drainage has.
Having defined the dimensionless groups, we also make all system variables dimensionless. We scale horizontal coordinates by the film half-length and vertical coordinates using the film's initial half-thickness. Moreover, we scale surfactant concentrations in the bulk and on the surface by surfactant concentrations in the bulk of the Plateau border and on the surface of the Plateau border, respectively. Velocities along the film are also nondimensionalized using a Marangoni velocity scale. Time is nondimensionalized using the ratio between the film half-length and the Marangoni velocity scale. The full details of the nondimensionalization can be found in supplementary Sec. S 5. In the dimensionless system, it turns out that the film half-length and the initial film half-thickness are set to unity, and surfactant concentrations in the bulk and on the Plateau border surface are also set to unity. Note that from now on, we only use dimensionless variables unless specified otherwise.

Velocity fields
A review of the derivation of dimensional governing equations is explained in supplementary Sec. S 3. However, here we present the equations in dimensionless form. The equations giving velocity fields in the x and z directions (the directions are as indicated in Fig. 1) are based on lubrication theory and turn out to be (Vitasari et al., 2013): where u; w are dimensionless velocities in x; z directions. Also, _ d is dimensionless rate of change of film half-thickness, d is dimensionless film half-thickness, and C is dimensionless surfactant surface concentration on the film. Note that w at z ¼ 0 is zero, while at z ¼ d, it equals _ d.
Since C might have an a priori unknown variation with x, Eqs.
(5) and (6) indicate that velocities are a priori unknown functions of position and time. Because we are dealing with surfactants that are being advected on the surface, we need to know about surface velocity. By setting z ¼ d in Eq. (5), we obtain an equation for surface velocity u s as follows: To proceed we also need to know the value of d at any given time t.
The Reynolds drainage formula for a rigid surface (Breward and Howell, 2002) gives in dimensionless form: where the respective solution for the film half-thickness is: Note the considerable simplifications (see Secs. S 2 and S 3 for details) that have gone into all these equations. Foam films are not only assumed to be thin (so lubrication theory applies), but in addition their thickness is treated as being spatially uniform, albeit varying in time. In reality foam film thickness can vary spatially as the foam film drains (Frankel and Mysels, 1962;Yeo et al., 2001;Joye et al., 1992). However describing that process involves the complication that the film surface becomes curved, and a normal stress boundary condition relating pressure to surface tension and curvature is then required at the surface. Here instead we simplify as per a previous study (Vitasari et al., 2013) by retaining only tangential surface stresses (due to the Marangoni effect associated with gradients of surfactant surface concentration). We also assume as per that previous study (Vitasari et al., 2013) that the film drains as if it were a uniform thickness squeeze film. Yet another simplification is that Eq. (5) and (6) are twodimensional equations (one dimension along the film, and one dimension across it). In reality flow is three-dimensional (two dimensions along the film and one dimension across it). The difficulty however is that the exact three-dimensional flow field that results depends on the film's shape, and specifically how many edges it has. Given that there will be considerable variation from film to film, not only in terms of film size, but also number of edges, we have for simplicity assumed a two-dimensional flow field here.
Having determined expressions for the velocity fields, we can now start to address the dimensionless mass transfer equations on and within the film surface.

Modelling of evolution of surfactant in a foam film
To investigate the evolution of surfactants, we use the dimensionless surfactant mass balance equations within the bulk and on the film surface. The dimensionless surfactant mass balance in the bulk follows a general mass transfer equation (Bird et al., 2007;Ubal et al., 2010;Cussler, 2009): @c @t þ r Á ðu cÞ ¼ 1 Pe where c is dimensionless surfactant concentration in the bulk, t is dimensionless time, u ¼ ðu; wÞ is dimensionless velocity in the bulk in x and z directions and D is the initial aspect ratio. A dimensionless diffusivity tensor (j) has appeared in the equation because of the nondimensionalization of the mass transfer equation using different characteristic lengths in x and z directions. Therefore, the coefficients of the diffusion terms in x and z directions are D=Pe and 1=ðPe DÞ, respectively.
As D is typically an extremely small number (i.e. the film is assumed to be very thin) and Pe is a relatively large number (see Table S 2), D=Pe is an extremely small number. This makes the diffusive flux in the x direction along the film negligible in all cases. Convection then dominates the mass transfer in that direction, and moreover since we have used a simplified flow field (as explained in Section 2.2), we have thereby simplified the convective mass flux. However, when Pe D is a reasonably small number, the diffusive flux in the z direction across the film becomes significant. Thus, small Pe D refers to the case in which there is a high diffusive flux across the thin film, even though the diffusive flux along it is negligible. This then is the case to be considered here, as will be discussed further in Section 2.4. The dimensionless mass transfer equation for surfactant in the bulk can now be simplified as: @c @t þ @ðucÞ @x þ @ðwcÞ @z Here c varies in the x direction, but only very weakly in the z direction due to Pe D being small. As well as accounting for surfactant in the bulk, we also must account for surfactant on the surface. Thus, we write general dimensionless surfactant mass balance on the surface as follows (Stone, 1990;Cantat and Dollet, 2012): where C is dimensionless surfactant surface concentration, u s is dimensionless velocity on the surface and S is the solubility parameter. Some previous works have ignored the term on the right hand side, which represents transport of surfactant from the surface into the bulk (Vitasari et al., 2013;Sonin et al., 1994;Jensen and Grotberg, 1992;Bruell and Granero-Belinchón, 2020;Gaver and Grotberg, 1990;Jensen, 1995). Here however that term must be retained. Nonetheless, one effect which is ignored in the present study is surfactant surface diffusion (Rio and Biance, 2014). The reason is that diffusion along the surface is generally much weaker than Marangoni flow on the surface. It is clear from the form of Eq. (12) that we cannot in general solve for C unless we also solve Eq. (11), and we know specifically how c varies with both x and z. Here however, as already alluded to earlier, we make a simplifying assumption namely the small Pe D limit, which allows us to progress the analysis even prior to solving for c. This is the subject of the next section.

Combining equations in the small Pe D limit
Here we model the evolution of surfactants on and within the foam film making use of the small Pe D assumption. To have a set of equations for the evolution of surfactant, regardless of whether in the bulk or on the surface, we combine mass transfer equations for the bulk and surface together. To do this, we integrate Eq. (11) over z from 0 to d, which results in: Notice here that in the final term on the left hand side, the term in c has been taken outside the integral sign. This is because as already mentioned, in the small Pe D limit, the value of c tends to vary only weakly across the film thickness. As Section 2.2 discussed, and as can be obtained from Eq. (6), w at z ¼ 0 is zero (because of the symmetry of the model), and w at z ¼ d is _ d. By combining two terms of the left-hand side of Eq. (13) and multiplying both sides of the equation by the S parameter, we obtain: The term on the right-hand side of the equation is the diffusive flux from the surface to the bulk, which is the negative of the diffusive flux from the bulk to the surface as appears in Eq. (12). The combination of Eqs. (12) and (14) leads to the following equation: From Eq. (15), we can see that the rate of surfactant change in the bulk and on the surface is related to the convection flows, both in the bulk and on the surface. However, we can only progress if we can identify a relationship connecting C and c. This can be done with an adsorption isotherm, which is the subject of the next section.

Adsorption isotherms
So far, we have presented the evolution of surfactant in the form of a single equation ( Eq. (15)). However, as this formula incorporates both C and c, we need an adsorption isotherm to link these two quantities together. Here we choose to express the isotherm in dimensionless form, which leads, as we will see shortly, to a rather simple relationship between these quantities C and c. The derivation of the original dimensional form and an explanation of how it is made dimensionless can be found in supplementary Secs. S 4 and S 5.
As Sec. S 4 explains, we start with a nonlinear adsorption isotherm called the Langmuir isotherm. This is a commonly used adsorption isotherm, which describes the relationship between surfactant surface and bulk concentrations (Chang and Franses, 1995;Hiemenz and Rajagopalan, 1997;Shchukin et al., 2001). The Langmuir isotherm is characterized by a so-called Langmuir parameter that quantifies the affinity for the adsorption and also by a maximum amount of surface concentration or surface excess. It is however possible to simplify the nonlinear isotherm to a linear one as we now go on to explain.

Approximations to the isotherm
Significant simplifications will arise in the governing equations if we manage to replace the Langmuir isotherm with a straight line.
To do this, we use two distinct approaches. First, we use a straight line to join the origin of a c; C plot (see Fig. 2)) to the conditions c Pb and C Pb at the Plateau border. This we call a global Henry isotherm approximation. In the other approach, we construct a tangent to the Langmuir adsorption isotherm, again at the Plateau border conditions, which we call a local Henry approximation. The detail of how we did this and a schematic representation can be found in supplementary Sec. S 4 and Fig. 2, respectively. Here as mentioned the approximate isotherms are to be used in a nondimensional form, such that c Pb and C Pb are in effect, both normalized to unity.
In this case, the global Henry isotherm in dimensionless form turns out to be: The local Henry isotherm is only slightly more complicated. We define a (dimensionless) local Henry constant (K HðlocÞ ) to be the slope of the (nondimensionalized) Langmuir adsorption isotherm at the point at which we construct the tangent, and it turns out to become: where K L is a dimensionless Langmuir parameter (see supplementary Sec. S 5 for how it is obtained). One of the features of a local Henry isotherm (evident in Fig. 2) is that even if we approach low concentrations, we have a minimum but still nonzero surface coverage. The minimum amount of surfactant surface concentration for a local Henry isotherm (C min ), in dimensionless form, becomes: Via Eqs. (17) and (18), the dimensionless local Henry isotherm becomes: Note that C min is just 1 À K HðlocÞ .
The dimensionless Langmuir parameter K L is what controls the value of K HðlocÞ , and as can be seen from Eqs. (17) and (18), when it is small (K L ( 1), the local Henry constant (K HðlocÞ ) becomes almost unity, and the minimum amount of surfactant surface concentration (C min ) becomes negligible. Hence, according to Eqs. (16) and (19), local and global Henry isotherms would become almost identical. On the other hand, if K L is large (K L ) 1), it follows that K HðlocÞ is small, whereas C min is close to unity. There is now quite some difference between the global and local Henry isotherms.

Considering dimensional variables for isotherms
The above discussion has been formulated, as we have said, in terms of dimensionless variables, with c Pb and C Pb in effect normalized to unity. However, it is also useful to consider what it means in terms of dimensional variables (such as Sec. S 4 considers).
Using a fractionation process with reflux tends to ensure that the surfactant concentration of the Plateau border is richer than in a system without reflux. In that case, and given that the purpose of fractionation is after all, to concentrate surfactant, the (dimensional) c Pb and C Pb are pushed towards the higher surfactant concentration part of the Langmuir adsorption isotherm plot, where surfactant surface concentration no longer changes significantly with changes in bulk concentration. In this situation, provided the concentration c in the foam film is not too far away from c Pb , then the tangent to the Langmuir isotherm at the Plateau border concentration or, in other words, the so-called local Henry isotherm describes the system's behaviour realistically (see e.g. Fig. S 1). On the other hand, a relatively dilute system could be considered, such that the surfactant concentration in the Plateau border is not exceedingly rich even in spite of reflux. In that case the (dimensional) c Pb and C Pb are pushed towards the lower surfactant concentration part of the Langmuir adsorption isotherm plot, where the slope of the isotherm is now rather larger. This incidentally tends to increase surfactant depletion length, which would decrease the value of S defined via Eq. (3). However, the relevant point here is that the Langmuir isotherm, global Henry isotherm and local Henry isotherm are all then rather similar, so the global Henry isotherm (which mathematically speaking is a little simpler than the others) would provide an acceptable approximation. To summarize we have defined so-called global and local Henry isotherms which, at least in certain parameter domains, are reasonable approximations to an overarching Langmuir isotherm. Although we could formulate a surfactant transport model utilising a Langmuir isotherm directly, we elect here to work with the global and local Henry isotherms instead. This is convenient as we will see, because it leads to a governing equation that is linear. Indeed in certain cases, i.e. without film drainage, we will see that it is even possible to obtain exact analytical solutions, in the form of series expansions. In any case, now that we have defined the link between the amount of surfactant on the surface and in the bulk and, we proceed to look at the total amount of surfactant (surface and bulk taken together) and how it evolves.
2.6. Determining total amount of surfactant at any film location When a fractionation system is used to recover surfactant, the distribution of surfactant between surface and bulk within the foam film is less significant than the total amount of recovered surfactant and how much the foamate is enriched. To begin to address this, we define a variable (C tot ), which is the total amount of surfactant present at any location in the film (or strictly speaking half the total amount at that location, because we consider just half of the thickness here). This can be expressed formally as: where C is dimensionless surfactant surface concentration, c is dimensionless bulk concentration, S is dimensionless solubility parameter and d is dimensionless film thickness. For the small Pe D limit, such that surfactants in the bulk are evenly distributed across the film, Eq. (20) can be simplified as: an approximation already used in Eq. (15). Note that only if we consider insoluble surfactant ðS ! 0Þ, does the value of C tot become almost equal to the amount of surfactant on the surface ðC tot % CÞ. For any other S > 0, the values of C tot and C differ. We can however, eliminate C and c from Eq. (15) in favour of C tot , making use of the earlier defined adsorption isotherms.
Moreover, once the value of C tot is known, it can be broken down into its separate components associated with C and c. To show how to do this, we use the local Henry isotherm (Eq. (19)) to derive the relevant equations. This is done because equations can be easily transformed to those for the global Henry isotherm if required, merely by setting K HðlocÞ ¼ 1 and C min ¼ 0. Via Eqs.
(21) and (19), we obtain: The combination of Eqs. (22), (23) and (15), together with in addition Eqs. (5) and (7) gives, after some algebra, the equation for the evolution of C tot . This equation becomes: Eq. (24) is a parabolic partial differential equation. Despite being linear, it is not always possible to solve analytically, given that coefficients depend on x and also implicitly on t, remembering here that d varies with time, at least when the film is draining. As a result, a method of solving this equation numerically has been applied, specifically a ''spectral method" (Brio et al., 2010). In this method, the solution of the equation is expressed as a finite expansion of some set of basis functions: we employed a Fourier series. Details of the numerical procedure are found in supplementary Sec. S 6.
When the film is not draining at all, so that d 1 and _ d 0, the situation is simpler still: each term in the Fourier series expansion then decouples and evolves independently of the others. A separation of variables solution then results, and so in effect, an exact analytical solution is obtained. This situation applies regardless of whether we employ a global or local isotherm, although the rate of evolution for each Fourier term is influenced by the choice of isotherm.
We also require initial and boundary conditions. We know that before the foam fractionation process starts, there is an initial total amount of surfactant on the film surface and in the film bulk which we denote C tot;0 and which is assumed to be uniform with x. Specifically C tot;0 is selected to reflect that there is less surfactant in the film initially compared to the corresponding amount in the Plateau border, remembering here that the Plateau border has been enriched by reflux. This difference between film and Plateau border provides (via the Marangoni effect) the main driving force for surfactant transport from the relatively surfactant rich Plateau border towards the relatively surfactant lean film. Note that given the value of C tot;0 , we can use Eqs. (22) and (23) to determine the corresponding initial surface concentration C 0 and initial bulk concentration c 0 . Alternatively we can use Eq. (21) to recover C tot;0 for given C 0 and c 0 . Remember that within all these equations, the value of d is unity initially (see Eq. (9)).
The boundary condition of the model at x ¼ 0 is dC tot =dxj x¼0 ¼ 0, which is due to the symmetry at the centre of the film. We also need a boundary condition at x ¼ 1 where the film meets the Plateau border. The challenges of imposing a boundary condition at this particular point have been discussed by Ruschak (1982Ruschak ( , 1985: even though flow in the film is a lubrication flow with near parallel streamlines, flow in the Plateau border meniscus itself is generally two-dimensional. Here in the interests of simplicity, we impose C tot j x¼1 ¼ C tot;Pb , where as already mentioned in Section 2.5.1, we assume throughout that the dimensionless c Pb and C Pb are fixed at unity. In other words, the Plateau border is treated as a reservoir of surfactant (Vitasari et al., 2013). The basis for this assumption however is that there might be rather more liquid in Plateau borders than in films, at least when films are thin: surfactant exchange from Plateau border to film is assumed to affect concentrations in the film, but has little impact on concentrations within the bulk of the border. A consequence of this though is that C tot;Pb , i.e. the total amount of surfactant at the location at which Plateau borders connect with films is actually weakly time-dependent if films gradually become thinner over time. Specifically it follows from Eq. (21) that C tot;Pb ¼ ð1 þ S dÞ where d is given via Eq. (9).

Determining overall amount of surfactant at any time
To assess the fractionation process performance, we need to calculate the overall amount of surfactant throughout the film (surface plus bulk) at each time. This is obtained by integrating C tot over the film half-length (which is unity in our dimensionless system). Thus In what follows, we call C ove the overall amount of surfactant, although strictly speaking it is the overall amount within a film half-length and half-thickness.
This quantity depends on time and hence upon the duration for which the film and Plateau border are in contact as reflux proceeds. This duration is dependent on process parameters, such as fractionation column length and velocity of bubbles through it (Shea et al., 2009). Thus, calculating C ove as a function of time can help to design or operate a more efficient fractionation column. Note that later on, we also compare the ratio between the overall amount of surfactant at each time and the overall amount of surfactant at the initial time, and we denote this C ove =C ove;0 . This measures the extent to which the surfactant recovered in the film is increased by contacting the film with a surfactant rich Plateau border. Note that because of the way we select the initial condition and also the way we normalize length in the dimensionless system, the value of C ove;0 is actually the same as C tot;0 .

Determining effective concentration at any time
Another parameter that helps us evaluate foam fractionation performance is a so-called surfactant effective concentration (c eff ). This quantity tells us how rich the foamate is at any instant. It can be expressed as the overall amount of surfactant on the surface and in the bulk at a specific instant, divided by the volume of the film. As we are dealing with a two-dimensional system with film half-length equal to unity, the volume can be just expressed in terms of the film half-thickness (d). This leads to Note that according to Eq. (21), the value of effective concentration c eff only matches the bulk concentration within the film in the limit when S d ) 1, i.e. for soluble surfactants and films that are not too thin. Otherwise c eff exceeds the bulk concentration.

Recovery and enrichment
Using the definitions in Sections 2.7 and 2.8 we can proceed to determine recovery and enrichment which are quantities often used to assess fractionation performance (Gerken et al., 2006;de Lucena et al., 1996;Exerowa et al., 2018;Stevenson, 2012). In this study, enrichment is defined as the ratio between the effective concentration of the foamate at each instant relative to the initial feed bulk concentration (c eff =c 0 ). Recovery on the other hand is taken to be measured just by the value of C ove , which specifically is the overall amount of surfactant recovered per film (or strictly speaking, per half-film length and half-film thickness). This differs from a conventional definition of recovery which would be the overall amount of surfactant recovered in the foamate relative to the overall amount fed to the fractionation column (Exerowa et al., 2018;Stevenson, 2012;Stevenson, 2014). Clearly, if we know the overall amount of surfactant recovered per film and also the number of foam films that enter the foamate, then we have a measure of the amount of surfactant recovered, which can be compared with the amount originally in the feed. Hence recovery per film C ove is proportional to the conventional definition of recovery, but it can be defined conveniently in the current model without needing to specify either the total number of foam films recovered or the overall quantity of surfactant in the feed. To summarize, in this work C ove represents recovery and c eff =c 0 represents enrichment.
As we will see later on, generally, there is a competition between recovery and enrichment (see Section 4.8), with high solubility parameter S favouring recovery and low solubility parameter S favouring enrichment.
Values of C tot ; C ove and c eff , as well as enrichment and recovery which derive from them, evolve with time and can approach a final state in the limit of very long times. This could either be a final steady state in the absence of film drainage, with an equilibrium between the Plateau border and film, or if there was film drainage, a quasisteady-state balance could be reached, as we discuss next.

Late time behaviour of the system under consideration
Recall that in Vitasari et al. (2013), the Marangoni flow is found to be the dominant contribution to the flow field at short times. However, at longer times, Marangoni tends to decay, such that a balance between Marangoni flow and film drainage flow then occurs, which leads the system into a quasisteady-state situation. Note that strictly speaking, this is indeed a quasisteady state rather than a true steady state, since the surfactant surface concentration continues to evolve, albeit slowly, even after Marangoni and film drainage flow balance. The difference in the present study is that surfactants are present in the bulk too. As a result, not only does the Marangoni-driven surface flow turn out to have slower impact on the surfactant surface concentration owing to surfactant escaping into the bulk (see Section 4.1 and later sections), but also film drainage removes liquid and hence surfactants from the bulk as drainage proceeds. Thus, although there is still an eventual balance between Marangoni flow and film drainage, this will happen at longer times compared to the insoluble surfactant case (Vitasari et al., 2013), and will also be further from equilibrating with the Plateau border (i.e. further from the state reached in the absence of film drainage). Note also that investigating the quasisteady state is primarily of interest in the case of the global Henry isotherm. This is because (as we will see later on, see Section 4.3 onwards) in the local Henry isotherm case, the impact of the Marangoni effect upon surfactant transport is slowed down even more significantly than in the global Henry isotherm case. Hence, for a local Henry isotherm, a balance between Marangoni flow and film drainage flow would likely only be realized at times much later than a typical residence time in a foam fractionation process. The following analysis focusses therefore on the global Henry case.
Note that thus far we have focussed primarily upon a partial differential equation for C tot (see Eq. (24)). Nonetheless using Eqs.
( 22) and (23) it is also possible to obtain analogous equations for C and c. If we are searching for a quasisteady state, it is easier to work in terms of C not C tot . This is because the boundary conditions for C are time-independent even with film drainage, whereas those for C tot are not. Specifically for C we find, still using the global Henry isotherm, Here Eq. (27) is the general unsteady equation for the evolution of C. It is necessary to explore whether or not this equation approaches a quasisteady state at a sufficiently long time. A description of the approximate analytical solution of Eq. (27) at late times can be found in supplementary Sec. S 9. In addition, further discussion and a comparison between analytical and numerical solutions are presented in Section 4.5. This completes setting up the model in the small Pe D limit. After deriving the equations for the evolution of C tot (or equivalently C or c), we need to use simulations to calculate C tot (or equivalently C or c) at every film position and every instant. These simulations are described next.

Simulation and benchmarking
Parameters corresponding to sodium dodecyl sulphate (SDS) have been used for simulation purposes within this study. SDS is one of the most studied anionic surfactants, and hence, relevant parameters are readily available in the literature, see Durand and Stone (2006) Tajima et al. (1970), Vollhardt and Emrich (2000), Wołowicz and Staszak (2020), Martínez-Balbuena et al. (2017), Nilsson (1957), Hines (1996), and Weißenborn and Braunschweig (2019). The origin of some of the important parameters used in our simulation and the related tables can be found in supplementary Sec. S 7.
Recall from Section 2.6 that a spectral method has been used to solve Eq. (24). Full details of the method are given in Sec. S 6. In our simulation, 38 Fourier terms have been used. Having generated the Fourier series, we can reconstruct the spatial variation of the solution from the Fourier series. To do this, we have evaluated at 500 spatial points throughout the film half-length. For the most part, 38 Fourier terms evaluated at 500 spatial points capture the surfactant distribution adequately, except for very slowly evolving systems at very early times (see e.g. Fig. 6 later on). Moreover, we have typically solved out to 20 dimensionless time units for the process duration divided into time steps of 0:001 dimensionless units. Time stepping is done via the fourth order Runge-Kutta method (Press et al., 1992). Note that based on the parameters we have chosen, each time unit is equivalent to 0:125 s (see Table S 1) thereby considering surfactant transport processes taking on the order of seconds. Note however that in certain cases (e.g. in Figs. 5 and 11 and also in Sections 4.5, 4.6 and 4.8), more than 20 time units have been considered to enable us to investigate the system's behaviour for somewhat longer times.
To benchmark our numerical method, we have compared our simulation with an alternative method called the ''material point" method, which had previously been used for the evolution of C in an insoluble case (Vitasari et al., 2013). As a result, this particular comparison has been made for the case in which the solubility parameter (S) is set to zero. These two numerical methods yielded almost the same results for C with a maximum relative difference of 3 Â 10 À5 . Details of the material point method can be found elsewhere (Vitasari et al., 2013).
We have carried out another benchmark using two analytical solutions for the case without film drainage: note that obtaining these analytical solutions relies on having a linear isotherm (either a global Henry isotherm or a local Henry isotherm, but not the Langmuir isotherm which is nonlinear). Neglecting film drainage then simplifies Eq. (24) and enables us to solve it analytically. The analytical solutions involve a complementary error function (erfc) approach or else (as mentioned in Section 2.6) a separation of variables solution approach, which results in a Fourier series solution. Both approaches are valid, but the complementary error function tends to be more convenient to use at early times, whereas the Fourier series is more convenient at later times. Details of these analytical solutions can be found in Vitasari et al. (2013). Further discussion of solutions in the absence of film drainage can also be found in Section 4, e.g. Sections 4.1 and 4.3.
Yet another benchmark used a case in which Marangoni flow has been switched off, but film drainage could still occur. This means that there is no mechanism driving surfactant transfer from the Plateau border toward the centre of the film. However, the transfer of surfactants from the film to the Plateau border still happens due to film drainage. This turns out to lead to an analytical solution, in which C tot is directly proportional to d, as can be easily verified from Eq. (24) upon dropping the Marangoni term. Comparing the spectral method with the analytical solution in the socalled ''no Marangoni" case has confirmed the consistency of the method once again.
Having thereby benchmarked the spectral method solution, we now consider the results that it predicts, including in cases which cannot be readily solved analytically.

Results and discussion
In this section, simulations for the evolution of the total amount of soluble surfactant (C tot ) have been carried out. As already alluded to, dimensional parameters for the system, assuming that the surfactant is sodium dodecyl sulphate (SDS), are themselves presented within supplementary Sec. S 7 and Table S 1. Dimensionless parameters are then reported in Table S 2. The solubility parameter (S, see Section 2.1) turns out to be 8:696, whereas the local Henry constant (K HðlocÞ , see Section 2.5.1) is 8:77 Â 10 À2 , and the film drainage velocity parameter (V R , see Section 2.1) turns out to be 0:0063. The global Henry constant is unity by definition. Other dimensionless parameters obtained for the system are also presented in Table S 2 including a check that the parameter Pe D really is smaller than unity.
To compare analogous conditions for different adsorption isotherms, we set the initial bulk concentration equal to half of the Plateau border's bulk concentration (c 0 ¼ 0:5) as the reference concentration. Therefore, C 0 is calculated to be 0:5 and 0:9561 for the global and the local Henry isotherms, respectively (Eqs. (16) and (19) and also Sec. S 8). In the local Henry isotherm case, note that C 0 is already quite close to unity, which has implications that we discuss later on. Regardless of whether we consider a global or local Henry isotherm case, once we know C 0 and c 0 , we now proceed to calculate the evolution of C tot . The evolution of C and c can also be calculated using Eqs. (23) and (22).
We begin this section by presenting the results for the case of the global Henry isotherm in the absence of film drainage in Section 4.1. Then, we discuss the global Henry isotherm case in the presence of film drainage in Section 4.2. After that, we consider the local Henry isotherm, first without film drainage and then with film drainage in Sections 4.3 and 4.4, respectively. An approximation to a quasisteady-state situation is also presented in Section 4.5. Details of how the quasisteady-state equations have been derived can be found in supplementary Sec. S 9. Then, in Section 4.6, we quantify the overall amount of surfactant in the film at any instant as well as the effective concentration, and compare the various soluble surfactant cases and the insoluble ones, with and without film drainage. Moreover, we investigate the effect of varying the solubility parameter upon the foam fractionation process: the results of this particular parametric study are found in Section 4.7. Finally, the performance of the fractionation column in terms of recovery and enrichment has been investigated and compared for different cases in Section 4.8.

Global Henry isotherm, no film drainage
We first neglect film drainage in the interests of simplicity, while also using the global Henry isotherm. It is reasonable to neglect film drainage as a first approximation, because the V R parameter is generally small. For the present case, within Eq.
Even though in this case Marangoni-driven transport is occurring along the film, it is accompanied by diffusion taking place across the film. Thus surfactant effectively escapes from the surface into the bulk of the film, and so does not accumulate quite so rapidly on the surface itself. Hence, compared to the insoluble surfactant case that has been presented by Vitasari et al. (2013), the rate of evolution of C tot is slowed down by a 1 þ S factor which is significantly greater than unity. This is as expected because the extent of the slow down depends upon how much surfactant enters the bulk, which in turn depends on solubility.
Note that the value of C tot at the edge of the film where it adjoins the Plateau border remains constant in this particular sys-tem without film drainage. The constant value is in fact just 1 þ S, and so is higher for a more soluble surfactant. The system also reaches a spatially uniform concentration C tot at very long times, equilibrating with the Plateau border. However, the time scale to reach this uniform concentration is much longer than in the insoluble case considered by Vitasari et al. (2013). When surfactant is soluble as in Fig. 3, even at 20 time units some nonuniformities are still evident, and equilibrium is still not reached.
This completes discussion of the case with a global Henry isotherm without film drainage. In what follows, we still consider the global Henry isotherm, but now, in the presence of the film drainage.

Global Henry isotherm, with film drainage
When film drainage is happening along with Marangoni flow, the surfactant transport process (see Fig. 4a) can be split into three stages of time. In the first stage at early time, there is a large concentration gradient near the boundary but no concentration gradient at the centre. Therefore, there is a strong Marangoni flow near the boundary which dominates the film drainage by a large amount. However, as there is still no Marangoni flow near the centre of the film, the effect of even comparatively weak film drainage can decrease slightly the value of C tot in this zone (Fig. 4a). It is noted further that, unlike the case considered in Section 4.1, there is now also continuous reduction in C tot at the Plateau border. This is just associated with the thinning process (see Eqs. (21) and (9) with C and c both unity at the Plateau border).
During the second stage occurring at intermediate times, the surfactant concentration gradient has now reached the centre of the film. The time scale needed for this to occur is however longer in the case of soluble surfactant than it is for insoluble surfactant (contrast the behaviour Figs. 4a and 4b): the reasons for this slower evolution in the soluble case have already been discussed in Section 4.1. During this second stage, the Marangoni flow driving surfactant onto the film is the dominant effect. This follows because there is still a significant difference between surface tension at the centre of the film and surface tension at the Plateau border, owing to the difference between the amount of surfactant in these zones. The film drainage, whilst present, is weaker than the Marangoni flow because of the slow film thinning rate.
In the third and final stage at later times, as the difference between concentrations along the surface becomes less and the Marangoni flow decays, film drainage is now comparable to the H. Rajabi and P. Grassia Chemical Engineering Science 265 (2023) 118171 Marangoni flow but acts in the opposite direction. This situation might lead to a quasisteady state, to be discussed further in Section 4.5. Particularly during this third stage, it is interesting to compare the soluble surfactant case in Fig. 4a with the insoluble case in Fig. 4b. When surfactant is insoluble, the amount at any spatial location always appears to increase with time as Fig. 4b shows.
As has already been mentioned above, in the soluble case however, this is not the case, as can be seen in Fig. 4a particularly when focussing on the region close to the Plateau border. At later times, because a significant amount of surfactant has transferred from the film surface to the film bulk, and because the film is itself draining, surfactant can actually be lost from the film.
Indeed the overall amount of surfactant in the film (C ove ) (see Section 4.6 for details) might start to decrease if the process continues for long enough. This would happen regardless of the amount of surfactant actually on the surface C, which always increases (see Fig. 5). Solubility might therefore actually have an adverse effect on overall surfactant recovered via foam fractionation if the film is draining and the process is left to run for too long in time, a point to which we return later on in Section 4.8.
Detail of what is happening to surfactant concentration C on the surface can be seen in Fig. 5. The advantage of plotting C instead of C tot is that, even with film drainage present, the value of C at the Plateau border does not depend on time. Hence an easier compar-ison between cases with and without film drainage can be performed. It is clear that in the presence of film drainage, the value of C is slightly lower than in the absence of film drainage. Moreover, in the case with film drainage, slight nonuniformities persist in C even at very late times: Fig. 5 extends up to 100 time units rather than just 20 units as in Fig. 4. Late time behaviour of the system is discussed further in Section 4.5.
This completes for now our analysis of the global Henry isotherm case. In what follows we switch to the local Henry isotherm, without film drainage in the first instance.

Local Henry isotherm, no film drainage
In the local Henry isotherm case (Eq. (19)), C min is not zero, and the initial amount of C (denoted C 0 ) for a given c 0 is rather higher than for the global Henry isotherm. As a result, C tot in the local Henry isotherm case starts initially from a slightly higher amount than for the global Henry isotherm. In addition, in our simulation, the local Henry constant (K HðlocÞ ) is determined to be 8:77 Â 10 À2 , a relatively small number (see Table S 2). This makes the surfactant concentration gradient on the surface and hence the Marangoni term much smaller than in the global Henry isotherm case. In the global Henry isotherm case, and without film drainage, the evolution has been slowed down by a 1 þ S factor because of the solubility. However, in the local Henry isotherm case, it is slowed down by a 1 þ S=K HðlocÞ factor (see Eq. (24)), which is a much more significant slow down. The issue is that there is now very little capacity to store additional surfactant on the surface, because C itself is already relatively high. Hence most of surfactant that is gained contributes to increasing the concentration c within the bulk of the film. Due to this very slow evolution, we see apparent oscillation at early times within Fig. 6. However, this is an artifact of using a limited number of Fourier components. Strictly speaking we need more Fourier components at those early times (Stroud and Booth, 2020;Boyd, 2000;Carslaw, 1925), but the oscillations soon decay away.
As can be seen in Fig. 6 (by contrast with Fig. 3), the local Henry isotherm starts from a slightly higher initial C tot than the global Henry isotherm does. This is owing to having a higher C 0 . However the main effect here is that the rate of increase of C tot is much slower than for the global Henry isotherm. Therefore, even after 20 time units, the Marangoni-driven surfactant transport has only just barely managed to reach the centre of the film.

Local Henry isotherm, with film drainage
Now, we consider the case of a local Henry isotherm in the presence of film drainage. In this case, again, having a small local Henry constant slows down the Marangoni flow. Apparent oscillations (which are numerical artifacts due to having a fixed number of Fourier terms) can again occur in this particular system, but to avoid them we simply started plotting from slightly later times. At the centre of the film, we can now see a reduction over time in C tot due to film drainage. This reduction is now much more significant than in Fig. 4a, which was obtained for a global Henry isotherm. The fact that a reduction might be seen suggests it may be important to manage carefully the foam film residence time in order to optimize the surfactant recovered, a point we return to discuss in Sections 4.6 and 4.8. Another possibility to consider to mitigate this is increasing the reflux flow through the system, which (see Section 2.1) thickens the Plateau borders and so can reduce the parameter V R and hence the film drainage rate.
The value of C tot is comprised (see Eq. (21)) of surfactant on the surface C plus surfactant in the bulk S d c, where recall S is given within Table S 2 and the evolution of d is given by Eq. (9). To understand these separate contributions to C tot , the evolution of C; c and d c are shown in Figs. 7b to 7d. It can be seen in Fig. 7b that the C values are always relatively high (i.e. close to unity), as the surface already has a significant amount of surfactant since the initial time.
Over time, there is a slight Marangoni-driven increase in C at positions neighbouring the Plateau border, remembering here that the value of C at the Plateau border itself is fixed. On the other hand, there is a slight decrease in C close to the centre of the film owing to film drainage.
Meanwhile, it can be seen in Fig. 7c, that changes in c are proportional to the changes in C. However, the range of c values encountered (from just less than 0:5 up to 1) is much larger than the range of C. However, we can see in Fig. 7d that it is d c which mirrors most closely the value of C tot . This reveals then that the H. Rajabi and P. Grassia Chemical Engineering Science 265 (2023) 118171 evolution of C tot , as obtained from Eq. (21), is significantly affected by both surfactant concentration in the bulk and by film drainage. To date we have analysed numerically the behaviour of surfactant on and within foam films up to some finite time. However it is also of interest to know how the system behaves in the limit of long times. Specifically, we want to obtain an approximate analytical solution for late times. This is the subject of the next section.

Quasisteady-state approximation for soluble surfactant
Here we discuss the expected long-time behaviour of the system, and compare it with the numerical results of the spectral method. The case without film drainage is simple because the film equilibrates with the Plateau border. Hence the case we consider here (as already alluded to in Section 2.10) is the one with film drainage.
The late-time behaviour for insoluble surfactant with film drainage has previously been discussed by Vitasari et al. (2013). In that study, it is explained that there is a quasisteady-state balance between Marangoni flow and film drainage terms at late times. A similar analysis can be carried out for a case with soluble surfactant (full details are in Sec. S 9). A complication is that (as we have already seen) soluble surfactant slows down the evolution of the Marangoni term, and hence delays the process of Marangoni and film drainage coming into balance. Under these circumstances (as already explained in Section 2.10), we choose to consider the global Henry isotherm case. In the local Henry isotherm case, as we have likewise explained, the evolution of the Marangoni term is slowed down even more. Marangoni and film drainage would then only balance after very long times indeed, which would likely exceed the residence time of a foam film within a foam fractionation process.
In a system containing soluble surfactants, when film drainage is active, the relevant parameter to consider when determining whether a quasisteady state is likely to occur turns out to be S V R . We know that V R is typically a small parameter while S is typically greater than unity. If S V R is a very large parameter, film drainage effectively dominates Marangoni flow (this case is similar to the no Marangoni case discussed in Section 3). The case of our interest is therefore, when S > 1, but S V R is still small compared to unity. The parameter values in Table S 2 indicate that this is indeed the case, so Marangoni flow and film drainage can then eventually come into balance. The derivation to obtain the approximate solution itself is presented, as we have said, in supplementary Sec. S 9.
Moreover, as Section 2.10 already explains, the approximate solution is expressed in terms of C rather than C tot , which is convenient because C has a boundary condition that is independent of time. This approximate solution is a series expansion in powers of the small parameter V R . When V R vanishes, the only steady state solution is to have C equal to unity, i.e. equilibration between the film and Plateau border. For small but nonzero V R , perturbations to C written as V R C 1 (first order correction) and V 2 R C 2 (second order correction) occur. Here C 1 and C 2 turn out to be functions of x and of d but not explicitly of time (although d itself depends on time, making the solution for C quasisteady). The value of C 2 also depends on the solubility parameter S, although it turns out that C 1 does not depend on S (see Sec. S 9 for details).
In Fig. 8, we can see a comparison between the results of the numerical solution and the approximate solution obtained in supplementary Sec. S 9. The data shown here correspond to 100 time units, for which a quasisteady state is likely to be a reasonable approximation, even in the case of a soluble surfactant. As can be seen, the values of C that have been predicted with the series expansion solutions are slightly larger than the values obtained numerically from the spectral method. However, including the second-order correction gives a better fit to the numerical solution result than including just the first-order correction does. The relative difference between C calculated from first-order correction and the one calculated using the numerical method at x ¼ 0, which has the greatest difference out of any x value, is equal to 5:47 Â 10 À4 , while the analogous difference for the second-order correction is 2:35 Â 10 À4 . It can also be seen that taking into account the solubility is required to have a more accurate approximation, as when the solubility parameter is ignored altogether, there is a much smaller difference that results between the firstorder and second-order corrections (see Sec. S 9.2 for details).

Quantifying overall amount of surfactant in the film
Although profiles of C tot ; C and/or c plotted against x as considered to date are relevant to study, what is of primary interest for the performance of a fractionation process is the overall amount of surfactant accounting for the entire film. In this section, overall amounts of surfactant, for different cases have been compared with each other. These cases consist of global and local Henry isotherms, with or without film drainage as well as soluble and insoluble surfactant cases.
The first comparison is for the overall amount of surfactant recovered per film (C ove ) calculated by Eq. (25). Another comparison is the ratio of the overall amount of surfactant in the film with respect to its initial amount (C ove =C ove;0 ): this is a measure of how much extra surfactant is recovered as a result of reflux. In addition, surfactant effective concentration (c eff ) has been calculated via Eq. (26) and compared in different cases: this measures how enriched the foam film is compared to an equivalent volume of bulk solution within the film itself. Fig. 9a shows C ove for different cases up to 100 time units. As is expected and can be seen from Fig. 9a, using a global Henry isotherm rather than a local Henry isotherm, leads to faster growth in the amount of surfactant at early times. However, in the global Henry isotherm case with film drainage, despite the comparatively fast initial increase, surfactant later decreases after reaching a maximum. A decrease at late times is also seen in the local Henry isotherm case, but there is now barely any growth in C ove at all before the decrease begins.
Another observation is that in the soluble cases, when film drainage is absent, the global and local Henry isotherm cases should eventually reach the same overall amount of surfactant. However, the local Henry isotherm case evolves slowly and so on the time scale of Fig. 9a still remains some way away from what is recov- Fig. 8. Profile of C calculated using a quasisteady-state approximation and numerically via a spectral method at 100 time units using a global Henry isotherm.
H. Rajabi and P. Grassia Chemical Engineering Science 265 (2023) 118171 ered the global Henry case. This implies that the surfactant recovered by the fractionation process would be sensitive to the residence time (or equivalently fractionation column length and/or bubble rise velocity). Note also that the lesser amount of surfactant predicted to be accumulated in the local Henry case has resulted when comparing the local and global cases at the same value of S. For a particular surfactant though, switching from a local to a global isotherm corresponds to lowering the target surfactant concentration (in dimensional variables) and, as has been explained already in Section 2.5.2, this could lead to an even lower S in the global Henry isotherm case. Finally, as is expected, the insoluble cases have the least overall amount of surfactant, there now being no surfactant in the bulk, while the differences between cases without and with film drainage cases are also very small. We use another comparison, the overall amount of surfactant at each time divided by the overall amount at the initial time (C ove =C ove;0 ), that can be seen in Fig. 9b. This quantity indicates, as we have said, the amount of surfactant recovered at any given time relative to the state of the films at initial time, and hence the impact that the reflux has had. As can be seen in Fig. 9b, although the insoluble surfactant cases have much smaller C ove at any given time compared to the soluble surfactant cases, their C ove =C ove;0 values grow more rapidly, with very little difference between cases without and with film drainage.
The next most rapidly growing C ove =C ove;0 occurs for the soluble surfactant case with a global Henry isotherm but without film drainage. The analogous value of C ove =C ove;0 for a global Henry isotherm with film drainage starts out similar but then peaks and decreases. Choosing a residence time close to the time of the peak would ensure that recovery is benefitting from reflux.
The slowest growing cases for C ove =C ove;0 are those for a local Henry isotherm, especially when film drainage is present, in which case C ove =C ove;0 barely grows at all. In the case without film drainage, C ove;0 is actually slightly higher for a local Henry isotherm than for a global Henry isotherm, which means that C ove =C ove;0 for the local isotherm at long times should always end up just slightly lower than for the global one. However, at 100 units of time the local Henry isotherm case is still quite some way from its final state.
The comparison of c eff as obtained via Eq. (26) in various cases can be seen in Fig. 9c, just considering cases with soluble surfactant. This gives a measure of how rich of the foamate is in surfactant, regardless of the actual amount of surfactant recovered. It is clear that longer residence times enrich the foamate, even if less surfactant is recovered. There is now relatively little difference between cases without and with film drainage, suggesting that even though drainage might be detrimental to the total amount of surfactant recovered, it is not detrimental to the concentration of what is recovered, and in fact sometimes it is beneficial.

Effect of changing solubility parameter
Recall that the main novelty of this study is the fact that we introduced the solubility parameter (S). We know (see Eq. (3)) that the solubility parameter is dependent on the initial film half-thickness and Plateau border's surface and bulk concentrations. Even though, in the particular system of interest here, we estimated the solubility parameter in Table S 2 to be 8:696, it is important to do a parametric study with respect to the solubility parameter. In Fig. 10, we used various solubility parameters to show the effect on the evolution of C ove =C ove;0 for both global and local Henry isotherms. Note that these comparisons have been made specifically in the case with film drainage, although unlike Fig. 9a we only considered times up to 20 time units.
By increasing the solubility parameter, the overall amount of surfactant in the film increases, but this effect is scaled out in Fig. 10 by comparing C ove =C ove;0 , rather than just C ove . It can be seen from Fig. 10a for the global Henry isotherm case, that systems with higher solubility parameter evolve more slowly, which is in line with expectations.
In Fig. 10b meanwhile, using the local Henry isotherm, Marangoni flow is slowed down very greatly, so the system is dominated by film drainage in the case when the solubility parameter is high. As a result, we lose surfactant from the film as the process evolves. Loss of surfactant from the film occurs sooner for higher values of the solubility parameter.
In Fig. 11, in addition to numerical data from the spectral method, we also plot the quasisteady prediction (see Sec. S 9) for C ove as a function of time for several S, with data extending now up to 100 time units. Eventually, the system does approach a quasisteady state. However, as shown in Fig. 11, in the cases with lar-ger solubility parameters, C ove calculated using the numerical solution only matches with the C ove calculated from the quasisteady prediction at rather late times.

Foam fractionation recovery and enrichment
Recovery and enrichment are important quantities that can be used to assess foam fractionation performance, see Stevenson et al. (2008Stevenson et al. ( , 2009Stevenson et al. ( , 2012Stevenson et al. ( , 2014Stevenson et al. ( , 2018. As Section 2.9 mentions, enrichment is the ratio between surfactant concentration in a foamate to the initial feed solution concentration and, in our model, is quantified as c eff =c 0 . Meanwhile C ove is used here to represent the recovery of the fractionation column at any instant. Again as Section 2.9 explains, strictly speaking this quantifies surfactant recovery per film, but it can be readily converted back to a more conventional definition of recovery once the number of films entering the foamate and the overall surfactant amount entering the feed are specified. Note that recovery and enrichment usually follow opposite trends (Saleh and Hossain, 2001;Gerken et al., 2006) and so, in a given system, to increase one we need to sacrifice the other. If we change the fractionation operation in some way however, we may be able to improve the recovery and enrichment performance. Changing residence time in the fractionation column and/or operating the fractionation with reflux as assumed here are possible ways of doing that. In fact, the recovery and enrichment at the initial instant in a column operated with reflux turn out to give a good representation of what the recovery and enrichment would be over a wide domain of times in an equivalent column without reflux. This follows because in the case without reflux, there is no Marangoni driving force tending to transfer surfactant from the Plateau border to the film. Any changes in recovery and enrichment over time then rely on film drainage which is comparatively weak. By the same argument, if operating with reflux is ever to be beneficial in recovery and enrichment terms, then it is necessary for residence time of foam films in a fractionation column to be sufficiently long for Marangoni flows to have taken effect.
In what follows, we present two sets of comparisons. First, we investigate the evolution of recovery and enrichment in the global and local Henry isotherm cases, both with film drainage. After that, we compare recovery and enrichment for the global and local Henry isotherms with and without film drainage, but just at selected times. In both sets of comparisons, different values of solubility S are considered across the domain 0:1 6 S 6 50.

Evolution of recovery and enrichment over time
As shown in Fig. 12a for the global Henry isotherm, there is a clear benefit in allowing the films to reside in the system for around 20 time units or more. The recovery vs enrichment curve has moved upward and to the right relative to the recovery vs enrichment in the initial state.
That said, the benefit of reflux is realised sooner in systems with small S, i.e. low solubility. In systems like that there is little advantage in extending the residence time much beyond 20 time units. At even longer times and for these low solubilities, we do see further slight improvements in enrichment due to film drainage, but the improvements are attained only slowly. On the other hand, for larger S values there is still benefit in extending residence time out to 50 or 100 time units. In systems with larger S, the Marangoni flow induced by the reflux itself impacts the system more slowly, and so may need longer residence times to take effect.
In the local Henry isotherm case (Fig. 12b), the situation is somewhat different. We now have more surfactant on the surface initially, which impacts the initial recovery vs enrichment curve, particularly at low solubilities for which systems are dominated by the surface rather than the bulk.
Indeed, across the full set of solubilities, the difference between the initial state and the state even at 100 time units is quite modest. We can divide the results plotted in Fig. 12b into three broad domains: low solubilities (on the bottom right of the figure), moderate solubilities (in the middle), and high solubilities (on the top left). In the low solubility domain (bottom right) there is little benefit of reflux at all. The system is dominated by the surfactant on the surface, but for the local Henry case, the surface concentration on the film already starts off very close to the surface concentration on the Plateau border, so barely changes over the course of time.  H. Rajabi and P. Grassia Chemical Engineering Science 265 (2023) 118171 In the moderate solubility domain (middle of Fig. 12b) we do see a gradual shift over time of the recovery vs enrichment curve, with longer times (even as long as 100 time units) being beneficial. Note that the equivalent solubilities in the global Henry isotherm case Fig. 12a had already converged after 20 time units, but in the local Henry case by contrast, the evolution is slowed down somewhat, hence the reason longer times are needed.
On the top left of Fig. 12b we see cases which even by 100 time units have barely shifted away from the initial state. This is because the effects of large solubility combined with a local Henry isotherm now slow down the Marangoni-driven evolution very significantly, meaning that reflux is yet to impact the system. To benefit from reflux at all, systems like this would need very long residence times indeed.

Comparison of recovery and enrichment of different systems at selected times
In this section, comparisons between different isotherms and with and without film drainage have been made. We examine the plots at just 20 and 100 time units.
We know that by construction there is no difference between cases with and without film drainage at the initial time. What Fig. 13a makes clear however is that there is still little difference between these cases even at 20 time units. Thus the main difference we see here at 20 time units is between the global and local Henry isotherm cases. This then manifests itself in the domain of larger solubilities (towards the top left of the figure). In the global Henry case, the film had already begun to acquire surfactant due to Marangoni-driven reflux (which is then beneficial for recovery and enrichment), but this has not yet happened in the local Henry isotherm case.
Meanwhile at 100 time units as can be seen in Fig. 13b, for less soluble surfactant (bottom right of the figure) there is little difference between the results for the two different types of isotherm. What we can see however is that cases with film drainage are being enriched slightly compared to cases without film drainage. Still at 100 time units but now for higher solubilities (top left of the figure), differences between the global and local Henry isotherm cases remain apparent. However, differences between cases with and without drainage are also clearly seen, with film drainage leading to less recovery. This is particularly evident in the local Henry isotherm case, for which loss of surfactant from the film via drainage has not been compensated by gain of surfactant via reflux-induced Marangoni flow.

Conclusion
Simulation of soluble surfactant transport on and within a foam film in the context of a foam fractionation column with reflux has been carried out. Reflux produces a Marangoni flow that drives surfactant onto the film. Despite incorporating tangential stresses that drive the Marangoni flows, the model remains nevertheless highly simplified: lubrication theory is used to determine flow fields, but films are assumed to remain uniform thickness, so film surface curvatures and normal stresses associated with them are not invoked. Likewise a two-dimensional flow field is assumed, recognising that the exact three-dimensional flow field that would otherwise result is sensitive to film size and also to number of edges a given film has. Variation from film to film is thereby expected.
The parameters used in the simulations were taken from relevant literature. As mentioned, this simulation assumes soluble surfactant, unlike previous work (Vitasari et al., 2013) which considered insoluble surfactant only. However, solubilities of different systems vary, and can be quantified by a solubility parameter (S). Moreover, a simplifying assumption has been applied in our modelling and simulations, namely, a small Pe D limit, in which a uniform surfactant concentration across the film thickness is obtained. This is what should happen in a sufficiently thin film when the diffusive transport across the film is fast compared to the Marangoni flow along it.
Our simulations consider two different adsorption isotherms: a global Henry isotherm and a local Henry isotherm. Both can be considered to arise from an overarching nonlinear isotherm (e.g. a Langmuir isotherm). The local Henry isotherm in particular is a new approach in which we use the local slope of the actual surfactant adsorption isotherm data and it gives a better approximation to the true isotherm in the higher range of concentrations. This is useful because in a typical foam fractionation process with reflux, the concentration in the Plateau border could well be situated at a point on Langmuir isotherm which is not too far from saturation, such that there is just a comparatively small change in surfactant surface concentration with respect to change in surfactant bulk concentration. As a result, a local Henry adsorption isotherm gives a better approximation in this domain of interest, while a global Henry isotherm still works reasonably well for significantly lower surfactant concentrations. Note also however that, for a given surfactant, decreasing surfactant concentration may also imply a decrease in the solubility parameter.
For each adsorption isotherm, two different cases are considered here, namely, no film drainage and with film drainage. Drainage, when included, is accounted for in a simplified fashion, i.e. still assuming the film remains uniform thickness but treating it now as a squeeze film. In the film drainage case, the amount of surfactant recovered by the film is less than in the case without film drainage, although the effective concentration of surfactant may be higher. The case without film drainage eventually equilibrates with the Plateau border, whereas the case with film drainage eventually approaches instead a quasisteady state. This quasisteadystate solution can be applied for later times when Marangoni flow and film drainage flow come into balance. The surfactant on and within the film then evolves only slowly with time due to a dependence of the quasisteady solution on the film thickness. Increasing the solubility parameter tends however to reduce the amount of surfactant on the surface in the quasisteady state, and moves the system further from equilibrium with the Plateau border.
There is however a question as to whether the quasisteady state is even reached in a typical residence time of a foam film in a fractionation column. Indeed, it has been found that the solubility reduces the impact of the Marangoni flow acting along the surface and so slows down the Marangoni-driven evolution. The reason is that the surfactant carried along the surface tends not to accumulate there but instead can escape into the bulk. This effect is particularly noticeable for the local Henry isotherm, as there is then little capacity for surfactant to accumulate on the surface, so a great deal necessarily escapes to the bulk. There also tend to be very low gradients in surfactant concentration along the surface in the local Henry isotherm case which slows Marangoni flows.
Moreover adding film drainage (which opposes Marangoni flow) slows the evolution even more compared to a case without drainage. Even though drainage is itself nominally weak, because it acts across the whole film thickness, its impact in mass transfer terms is not necessarily less than that of the Marangoni flow.
Higher solubilities give at any specified time, less growth in the overall amount of surfactant relative to its initial amount. This is due to the previously mentioned fact that higher solubilities slow down the Marangoni flow. Moreover, as film drainage is now in relative terms more important, there will be a bigger deviation at later times from conditions applicable in the Plateau border. The film drainage might even cause the reduction of the total amount of surfactant recovered after a certain time, particularly in cases with significantly slowed down Marangoni flows. This then impacts on the amount of surfactant recovered. However, despite the decreasing recovery, the enrichment always increases because both film drainage and Marangoni flow contribute to having richer films.
In summary, the results of the study can help to understand the evolution of surfactants on and within foam films during foam fractionation with reflux. It is seen that reflux is beneficial for fractionation even of soluble surfactants. However, systems exhibiting a global Henry isotherm (or an isotherm approximating to one) benefit more from reflux than systems with a local Henry isotherm (or likewise an isotherm approximating to one). Indeed the local Henry isotherm case would only benefit from reflux if residence times are very long. This knowledge can help to improve design and operation of fractionation columns.
Although an often overlooked element, i.e. the solubility parameter S, has been included in the model considered here, other improvements to the model are nonetheless still possible. For instance, the small Pe D assumption employed here is a reasonable assumption but would only be true in a sufficiently thin film. The value of the diffusivity coefficient is moreover another factor that might adversely affect the applicability of the small Pe D assumption. Generally, bulkier molecules tend to have smaller diffusivity coefficients. Therefore, Péclet number and hence Pe D for those bigger molecules (e.g. fractionation of a protein) can be much greater than for a comparatively small molecule such as SDS as has been considered here. In addition, by introducing the local Henry isotherm, we have tried to use the simplicity of a linear isotherm and simultaneously improve its accuracy over the global Henry isotherm for a higher range of concentrations. Nonetheless taking in account the full nonlinearity of the adsorption isotherm would be another means to improve this model. Additional improvements would involve relaxing some of the simplifying assumptions already mentioned above, e.g. allowing film thicknesses to be nonuniform rather than uniform, and allowing flow fields to be threedimensional rather than two-dimensional.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.