Transport of Convected Soluble Surfactants on and within the Foam Film in the Context of a Foam Fractionation Process

This study models convective transport of soluble surfactant in a foam fractionation system with reflux. Owing to reflux, initial surfactant concentration in films is lower than in Plateau borders. Marangoni flows and film drainage flows arise convecting surfactant both on the film surface and in the bulk. An interface is set up within the film bulk called a separatrix: this divides two regions of uniform surfactant concentration, one with concentration equal to that of the initial film and one with concentration equal to that in the Plateau border. The evolution of the separatrix is tracked to determine surfactant recovery and enrichment for the film. Surfactant lean films benefit most from reflux. However for films that are initially comparatively surfactant rich, recovery might actually decrease at long times owing to film drainage. Nonetheless surfactant lean films and those containing surfactant with only moderate solubility benefit from reflux even despite film drainage.


Introduction
Foam separation techniques have been identified as alternatives to more conventional separation processes such as ion exchange and ultrafiltration (Wong et al., 2001), particularly because of their efficiencies in dealing with dilute aqueous systems (King, 1980). Due to these advantages, these methods have recently been finding their place in various sectors such as pharmaceutical, environmental-related and biochemical industries (Boyles and Lincoln, 1958;Gehle and Schügerl, 1984;Grieves and Wang, 1967a,b;Husband et al., 1994;Lee and Ryu, 1979;Shao et al., 2020;Vitasari et al., 2013a;Xanthopoulos and Binnemans, 2021). One of the foam separation methods which is the subject of the present study is foam fractionation. Foam fractionation is a physicochemical process in which surface-active chemicals are separated from an aqueous solution by adsorption on bubbles rising in a column (Lemlich, 1972). As there are no solvents other than water existing during this process, it can be considered as a "green" process in sustainability terms (Burghoff, 2012). Some of the advantages of foam fractionation are low capital cost, low energy requirement and subsequently low operational costs (Wong et al., 2001). As a result, there have been various studies carried out in this field (Brown et al., 1990;Buckley et al., 2022;Burghoff, 2012;Du et al., 2000Du et al., , 2002Grassia, 2023;Hutzler et al., 2013;Keshavarzi et al., 2022;Liu et al., 2017;Mukhopadhyay et al., 2010;Shedlovsky, 1948;Tobin et al., 2014).
Some of these studies show a potential benefit of foam fractionation by returning part of the so called foamate back into the fractionation column, known as reflux (Lemlich and Lavi, 1961). Reflux, effectively puts into contact a rising stream of leaner bottom solution and a falling stream of enriched collapsed foamate (Martin et al., 2010). Rich liquid is then travelling through a network of Plateau borders contacting the foam films. Thus, there is an opportunity for further enrichment of the foamate (Brunner and Lemlich, 1963;de Lucena et al., 1996;Lemlich and Lavi, 1961;Martin et al., 2010;Rajabi and Grassia, 2023;Rubin and Melech, 1972;Stevenson and Jameson, 2007;Stevenson et al., 2008;Vitasari et al., 2013b).
Plateau borders referred to above specifically are liquid channels embedded within the foam where three films meet (Grassia et al., 2001;Weaire and Hutzler, 1999), and as we have said, they form a network. Due to often having a higher liquid volume compared to the liquid in the films adjacent to them, Plateau borders may to an extent be considered as surfactant reservoirs (Vitasari et al., 2013b). Plateau borders can therefore in principle transfer significant surfactant to adjacent films as reflux proceeds. However to achieve this surfactant transfer, flow must occur from Plateau border to film. For the purpose of this work, we use the term "surfactant" to encompass any surface-active molecule including big, bulky molecules like proteins which are often targeted for separation by fractionation (Brown et al., 1990;Du et al., 2000Du et al., , 2002Mukhopadhyay et al., 2010). These then are the molecules that, for reflux to be effective, must flow from Plateau border to film.
Although, in line with previous work (Rajabi and Grassia, 2023;Vitasari et al., 2013b), we focus the discussion here on foam fractionation with reflux, the results are also relevant to foam fractionation operated in another mode, namely stripping mode (Datta et al., 2015;Kamalanathan and Martin, 2 2016; Rubin and Melech, 1972). In stripping mode a feed is provided to a fractionation column and flows down through Plateau borders, contacting foam films as it flows. As with reflux, surfactant is again transferred from Plateau borders to films. However the objective now is not so much to enrich the films, but rather to remove surfactant from the liquid in the Plateau borders. This mode of operation would for instance be relevant for removing a surface active contaminant from a wastewater stream. For the most part in what follows, for simplicity we discuss reflux, even though stripping is also relevant.
Despite the potential advantages of foam fractionation over other separation methods, in view of the complicated flows that arise in the presence of reflux, the process of foam fractionation with reflux requires more research. In particular, it is useful to have a modelling study which can predict the extent to which reflux permits foam films to become enriched in surfactant, under different sets of conditions encountered during fractionation. This can subsequently help us to design and operate a more efficient fractionation column in the future. Specifically the model to be used in the present work has been built upon two previous modelling studies on film scale surfactant transport during foam fractionation with reflux (Rajabi and Grassia, 2023;Vitasari et al., 2013b), and in what follows we review them.
The main effect that was included in those studies was Marangoni flow taking surfactant rich material from Plateau border into the film. Moreover, film drainage which thins the film and which causes a flow towards the Plateau border, opposes Marangoni flow on the surface. In Vitasari et al. (2013b), the authors worked out the evolution of an insoluble surfactant over just the surface of a foam film. They also discussed a so called quasi-steady state for cases in the presence of film drainage, in which a balance on the film surface eventually happens between Marangoni flow and film drainage. However, the fact that surfactants have some level of solubility was neglected. Solubility must however affect the transport behaviour, which thereby affects the foam fractionation process.
Later, Rajabi and Grassia (2023) considered the solubility of surfactants within a foam film in addition to their presence on the surface. Marangoni-driven and film drainage-driven surfactant transport again occur, but transport now occurs not just on the foam film surface, but in the bulk of the foam film as well. However the study of Rajabi and Grassia (2023) focussed on a particular limit in which surfactant also diffuses rapidly across the foam film. This limit could be quantified in terms of dimensionless groups corresponding to a small Pe ∆ number, where Pe denotes Péclet number (measuring the ratio between convective and diffusive transport) and ∆ is the film thickness to film length aspect ratio. It was identified by Rajabi and Grassia (2023) that this particular limit would be a reasonable approximation for smaller surfactant molecules (which tend to have comparatively high diffusivity coefficients) being transported across particularly thin films. What was found by Rajabi and Grassia (2023) is that the impact of Marangoni flow is slowed down due to the solubility, compared to the case considering insoluble surfactants. The reason for this slow down was found to be surfactant escaping into the bulk of the film, once the Marangoni flow transported it from the Plateau border onto the film. It was also confirmed that the quasi-steady state condition introduced by Vitasari et al. (2013b) can also occur, after a sufficiently long time: Marangoni-driven and film drainage-driven transport are then effectively in balance.  However, there could be an opposite limit in which Pe ∆ is relatively large. This limit is relevant in cases with relatively small diffusivity coefficient, e.g. for bigger molecules such as proteins . It could also be relevant for foam films shortly after they are formed, such that the film is still in the process of draining, and hence could be quite some way from a final thickness. This limit is the subject of and the novel contribution of the present study (see Fig. 1). Diffusion across the film is now relatively slow and, as a consequence, transport from the film surface to the film bulk is relatively slow. Hence, what happens on the surface turns out to be identical to what has been discussed by Vitasari et al. (2013b) (albeit distinct from the system considered by Rajabi and Grassia (2023) which has small rather than large Pe ∆). Flow is affected by competition between Marangoni-driven and film drainage-driven transport on the surface. However the Marangoni and film drainage flow fields are not just confined to the surface but exist in the bulk as well. Since, similar to Rajabi and Grassia (2023), surfactants are now treated as soluble and present in the bulk, the aforementioned flow fields necessarily convect surfactant in the bulk. In the large Pe ∆ limit considered here, convection (not diffusion) is the dominant transport mechanism in the bulk. Indeed, as we will see, convection plays an important role in transferring surfactants from the bulk of the Plateau border into the bulk of the film.
A point to emphasise is that, since the flow on the surface here is identical to what was found in the prior work by Vitasari et al. (2013b), the flow field in the bulk also turns out to be the same. In effect therefore the results to be obtained here could have been obtained by post-processing the results of Vitasari et al. (2013b). However such post-processing was never attempted by Vitasari et al. (2013b): for an insoluble surfactant as considered in that work, any motion in the bulk is irrelevant to surfactant transport. The novelty of the present work is therefore to consider soluble surfactant, for which transport in the bulk is certainly relevant, surfactant in the bulk now being passively convected.
In the model considered here, we start with a situation in which the Plateau border is surfactant rich and the film is surfactant lean (Fig. 1a). Marangoni stresses offset by film drainage set up a convective flow field. Convection then carries surfactant rich material from Plateau border to film and carries surfactant lean material from film to Plateau border (Fig. 1b). The net effect is that the film is enriched. This then is the process we wish to model.
Note that what we are trying to describe here (analogously to the work of Rajabi and Grassia (2023); Vitasari et al. (2013b)) is specifically the microscale process of surfactant transfer between an individual Plateau border and an adjacent foam film. We are not endeavouring at this stage to incorporate this microscale behaviour into a fractionation model at the entire process scale. In a typical fractionation column there will be a multitude of foam films and a multitude of Plateau borders, and the surfactant content of each depends on where in the fractionation column they are located, and also on the time elapsed since the fractionation process started. Sufficiently long after start up, a continuous fractionation system with reflux should settle into a steady state operation. However the time scale for that to happen presumably relies on residence times of flowing foam films passing up through the entire column and likewise residence times of flowing liquid (within Plateau borders) passing down through the column. This is not however what is focussed upon here.
The intention here is to keep the models as simple as possible, while still endeavouring to capture the main transport processes that are expected to occur. The models used by Rajabi and Grassia (2023); Vitasari et al. (2013b) were likewise greatly simplified, and analogous simplifications will be employed here (see also Sec. 2.3 for further discussion of some of the simplifications employed).
Two-dimensional rather than fully three-dimensional flow fields will for instance be considered. Given however that films in a real foam have different shapes and sizes, capturing the full threedimensional geometry of the flow on each and every film would be challenging in any case. Films are also to be treated as having a spatially uniform thickness, albeit with that thickness changing over time. There are of course models in the literature that study the fluid mechanics of a draining foam film in a much more sophisticated fashion accounting for non-uniformities in thickness (see e.g. Frankel and Mysels (1962); Joye et al. (1992); Yeo et al. (2001)): film surfaces are no longer flat leading also to pressure jumps across them. The film shape and the surfactant mass transport must then be determined together. These sorts of complications have been neglected in previously mentioned work on foam fractionation (Rajabi and Grassia, 2023;Vitasari et al., 2013b) and will be neglected again here also. Indeed we assume the film geometry and how it evolves is known, and focus just on surfactant mass transport. As we have alluded to though, the models to be used here still capture key physics, such as Marangoni flow, drainage flow and convective surfactant transport.
One aspect that will however differ here from previous work are the physicochemical parameters that we assume. Previous work (Rajabi and Grassia, 2023;Vitasari et al., 2013b) utilised parameters relevant to the common surfactant SDS. For the simulations here however, we use parameters relevant to the protein beta-lactoglobulin (β-LG) (Xiong and Bansal, 2022). As a much bigger 5 molecule than SDS, this has lower diffusivity , and as a result can readily meet criteria for having large Pe ∆ which is the basis of our model. Physicochemical data for β-LG are readily available, since it is a widely studied protein in foam and interface science (Atkinson et al., 1995;Baeza et al., 2005;Bouyer et al., 2013;Fainerman et al., 2020;Gochev et al., 2013Gochev et al., , 2021Gochev et al., , 2019Krägel et al., 2003;Lech, 2016;Lexis and Willenbacher, 2014;Miller et al., 2004;Pradines et al., 2009) and it has also been used in the context of foam fractionation (Cheang and Zydney, 2003;Shea et al., 2009). As has been mentioned though, we will often use the generic term "surfactant" to keep the discussion more general: the model we present requires large Pe ∆ but is not tied to any specific material, provided parameter values are available.
This study is laid out as follows. Sec. 2 outlines the mathematical theory used in the study of convected soluble surfactant transport, which is based on adapting the works of Vitasari et al. (2013b) and Rajabi and Grassia (2023) to this new convection-dominated system. In Sec. 2 equations are mainly presented in dimensionless form, but the nondimensionalisation process itself is presented in Sec. S 1 of the supplementary material. Technical details of numerical algorithms employed are also relegated to supplementary material. Algorithms already used by Rajabi and Grassia (2023) can be adapted supplemented with some additional features, so any discussion of algorithms in the supplementary material is focussed on those additional features. In particular the challenge of carrying out the calculations at early times is discussed Sec. S 2. Other than that, the numerical approach is similar to what Vitasari et al. (2013b) and Rajabi and Grassia (2023) have already done. However here algorithms take account also of mass transfer by convection in the bulk of the foam film, and a discussion of that can be found in Sec. S 3 of the supplementary material: convection in the bulk did not need to be addressed in the work of Vitasari et al. (2013b), but it certainly must be considered here, so the discussion of Sec. S 3 covers that. In addition, Sec. S 4 deals with selecting simulation parameters and benchmarking. Returning to the main text, results and discussion are presented in Sec. 3, with some supplementary results in Sec. S 5. Finally Sec. 4 deals with conclusions of the study.

Mathematical model for convected soluble surfactant transport
This study models convected soluble surfactant transport on and within a foam film. We first present in Sec. 2.1 essential dimensionless groups used to carry out the study. Derivation of the velocity fields in the bulk and on the surface of the foam film are discussed elsewhere (Vitasari et al., 2013b). However, due to their importance within the current study, they are also mentioned here in brief (Sec. 2.2). Mass transfer equations are presented in Sec. 2.3. As will be discussed in detail later, due to the Marangoni-driven and film drainage-driven convective flows in the bulk of the film, two regions form within the bulk. One region remains at the initial bulk surfactant concentration and the other region has a higher bulk concentration corresponding to the concentration found in the bulk of the Plateau border. The novel contribution of the present work is to model, simulate and analyse the evolution of these separate regions in the bulk and the subsequent effect on the fractionation process performance. We assess the size of these two regions by examining a cross 6 section of the foam film, looking at the total area of the cross section of the film and also the areas corresponding to each of these regions. The details of computing these two regions and the surfactant they contain can be found in Secs. 2.4 and 2.5. Finally Sec. 2.6 deals with overall measures of fractionation performance, namely recovery and enrichment.

Dimensionless groups
Definitions of dimensionless groups that have been used in this study are identical to the ones used by Rajabi and Grassia (2023). However here we are looking at a different parameter regime. Hence we present the dimensionless groups in brief.
The first dimensionless group is Péclet number. It is obtained (Leal, 2007) based on balancing surfactant transport by Marangoni effects along the film in the x-direction, and diffusive transport of surfactant across the film, in the z-direction (the directions are as indicated in Fig. 1), and can be expressed as follows: where G is Gibbs parameter (that measures sensitivity of surface tension to surface concentration) (Lucassen-Reynders et al., 2001, δ 0 is initial film half-thickness, µ is liquid viscosity, L is film half-length and D is diffusion coefficient (in this case for β-LG). In fractionation applications, Péclet number is typically a relatively large number, especially for proteins that due to their bigger molecule sizes have smaller diffusion coefficients (Tang et al., 2022) (see Table S 3 in the supplementary material for a typical value of Pe). Meanwhile, ∆ is the initial aspect ratio between film half-thickness and film half-length, and is defined as: where δ 0 is initial film half-thickness, and L is film half-length. In a typical foam film ∆ is a relatively small parameter (see Table S 3). However here, as can be seen in Sec. S 4.1 (and in particular in Table S 3), in the particular limit of interest, the product of Pe and ∆ remains a relatively large number which is in line with our assumptions.
Another parameter relevant to this study is the solubility parameter S. Solubility parameter describes the typical amount of dissolved surfactant relative to the amount of surfactant on the surface, and is defined as below: where δ 0 is initial film half-thickness, Γ Pb is surfactant surface concentration at the Plateau border, and c Pb is surfactant bulk concentration at the Plateau border: both Γ Pb and c Pb are dimensional quantities here, although later on we will also define dimensionless analogues of them. Note that for a given Γ Pb an insoluble or almost insoluble surfactant will have a very small c Pb and hence a very small S, but a more soluble surfactant will have a larger c Pb and hence a larger S value.
Formally Γ Pb /c Pb is a depletion length (extent of a bulk region containing an equivalent amount of surfactant as the surface itself), and S is then the ratio of the actual geometric extent to that depletion length. As Table S 3 makes clear, in systems of interest S turns out to be a dimensionless parameter on the order of magnitude of unity, and we will allow it to vary during the course of this study.
Finally dimensionless film drainage velocity parameter (V R ) is the ratio between velocity of film drainage under the action of capillary suction and typical velocity of Marangoni convection at the start of the process, which can be expressed as follows: where δ 0 is the initial film half-thickness, γ Pb is surface tension at the Plateau border, a is the Plateau border's radius of curvature and G is Gibbs elasticity. As Table S 3 makes clear, V R is typically a small parameter. At least early on in the process then, film drainage is a weaker effect than Marangoni flow. Moreover since the drainage velocity parameter is low (in other words, film drainage flow is slow), effects of drainage require quite some time before impacting the system.
To summarise, Table S 3 gives the typical values of dimensionless groups, and is based on dimensional parameter values obtained from literature (see Tables S 1 and S 2).
In line with what has been done in Vitasari et al. (2013b) and Rajabi and Grassia (2023), in addition to defining the above dimensionless groups, we also make all the system's variables dimensionless. The process of nondimensionalisation is again similar to Vitasari et al. (2013b) and Rajabi and Grassia (2023) and can be found in Sec. S 1. Briefly, we scale horizontal coordinates by the film half-length and vertical coordinates using the film's initial half-thickness. In addition, the crosssectional areas within the bulk of the foam film are nondimensionalised with respect to the initial film area. Meanwhile, we scale surfactant concentrations in the bulk and on the surface of the film by surfactant concentrations in the bulk and on the surface of the Plateau border, respectively. Velocities along the film have also been nondimensionalised using the Marangoni velocity scale, but transverse velocities have an additional factor of aspect ratio included in the scaling. Time is nondimensionalised using the ratio between the film half-length and the Marangoni velocity scale. Note from Table S 2 that the characteristic time scale is actually rather short, significantly shorter in fact than was the case in Rajabi and Grassia (2023), owing to focussing on not quite so thin films in the present work which admit higher velocities along them. Certainly typical film residence times in a foam fractionation column are likely to be many dimensionless time units.
Note that from now on, we only use dimensionless variables unless specified otherwise.

Velocity fields
To work out the velocity fields in the bulk and on the surface, we use the same approach taken by Vitasari et al. (2013b). The physical mechanisms that drive the flow fields are Marangoni stresses 8 and film drainage. Details of these mechanisms are discussed elsewhere (Rajabi and Grassia, 2023).
To derive the velocity fields, a lubrication approximation has been used which leads to the dimensionless velocity fields in the bulk as follows: where u and w are dimensionless vertical and horizontal velocity components in the bulk and Γ is dimensionless surfactant surface concentration. Here also δ is dimensionless film half-thickness andδ is dimensionless rate of change of δ with time. By the same token and using Eq. (5) with z = δ, the velocity field u s on the surface becomes: Note that in the present model, even though we are considering a soluble surfactant, the velocity fields here are unchanged from those considered by Vitasari et al. (2013b) for an insoluble surfactant. This is because (unlike the work of Rajabi and Grassia (2023) which allowed surfactant to diffuse readily between bulk and surface), here it turns out that surfactant fails to diffuse off or onto the surface on time scales of interest. Thus as far as setting up the velocity field is concerned, the surfactant might as well be insoluble. What is different from Vitasari et al. (2013b) is that this same velocity field, established as a result of conditions on the surface, now transports surfactant both on the surface and in the bulk. Hence, we move on to the mass transfer equation.

Mass transport equation in a foam film
The general dimensionless mass transfer equation for surfactant in the bulk can be expressed as follows (Bird et al., 1960): where, c is surfactant concentration in the bulk (which recall is made dimensionless here with respect to the analogous concentration in the Plateau border), Pe is Péclet number and ∆ is initial aspect ratio. In the limit of interest of the present study, Pe ∆ (see Table S 3) is large and Pe /∆ is extremely large. Hence diffusion terms in both x and z directions are negligible. Moreover, using the continuity equation for an incompressible liquid, Eq. (8) turns out to give Dc/Dt = 0, which means that following an element of fluid, there will be no change in the concentration with time.
The general mass transfer equation on the surface is as follows (Rajabi and Grassia, 2023): where Γ is surfactant concentration on the film surface (again made dimensionless using the analogous Plateau border surface concentration) and S is the solubility parameter defined earlier. The term on the right-hand side of the above equation is the diffusive flux from the bulk to the surface (in the z direction). It turns out to be negligible due to the large Pe ∆ assumption (with S being order unity here). This then confirms that, at least on the time scale of interest for Marangoni flow and film drainage, there is not any diffusive transport from surface to bulk. As Eq. (9) shows, the larger the value of S, the more likely it is that bulk-to-surface transport becomes relevant, but here Pe ∆ is much larger than S, so convection along the surface is much faster than any bulk-to-surface transport (which ultimately requires diffusion to be active).
As already alluded to, as far as mass transport on the surface is concerned, the system is then equivalent to the insoluble case which has already been considered by Vitasari et al. (2013b). Combination of Eq. (9) and Eq. (7) and using the mentioned assumptions in Vitasari et al. (2013b) leads to the following mass transfer equation on the surface: Eq. (10) has initial and boundary conditions defined as follows: where in particular Γ 0 < 1 so that the film starts off leaner than the Plateau border. The value of δ evolves according toδ = −V R δ 3 (Rajabi and Grassia, 2023) and hence where recall V R is a relatively small parameter, so δ evolves comparatively slowly.
As can be seen from Eq. (10), the evolution of surfactant on the surface is due to a competition between Marangoni flow and film drainage. Eq. (10) is a parabolic partial differential equation. A method of solving this equation numerically has been applied, specifically a "spectral method" (Canuto et al., 2006) (details can be found in Rajabi and Grassia (2023)).
As was mentioned previously, this study uses a similar set of simplifying assumptions used previously by Vitasari et al. (2013b) and Rajabi and Grassia (2023) (two-dimensional system, lubrication approximation, film surface remains flat, etc.).
One significant simplification that we highlight, is the assumed boundary condition in Eq. (11) that the Plateau border surface remains at (dimensionless) surface concentration unity. To justify this, it is important to recall a physical picture of how foam fractionation with reflux operates. Foam films rise up through the fractionation column, whilst simultaneously reflux liquid drains down through a network of Plateau borders, and we are looking at the mass exchange process between the two. For sufficient reflux flow, it can be the case that there is more liquid in the Plateau borders than in the films (Vitasari et al., 2013b). Similarly there can be more liquid flux (and also more surfactant flux) in the Plateau borders than in the films (Grassia, 2023). As a result, the relative change in surfactant content in the comparatively surfactant lean films should be greater than the relative change in surfactant content in the surfactant rich Plateau borders. Even though the Plateau border surfaces do lose surfactant to the films, the surfactant remaining on them is carried down under gravity, and so is replaced by additional surfactant arriving from higher up. In effect therefore Plateau borders are approximated here as being surfactant reservoirs. Ultimately though what the boundary condition in Eq. (11) attempts to capture is the notion that Plateau borders are richer in surfactant than films.
Here of course since surfactant is soluble, we do not consider just surfactant transport on surfaces, but also surfactant transport in the bulk. That said, even though the surfactant is soluble here, we reiterate that there is no transfer between bulk and surface on the time scale of interest due to the large Pe ∆ assumption. Convection between the bulk of the film and the bulk of the Plateau border is still permitted, but to the extent that the Plateau border is treated as a surfactant reservoir, what we must focus upon here is surfactant convection in the bulk of the film. Treating a Plateau border as a reservoir, as is done here, is arguably more relevant for foam fractionation with reflux than for foam fractionation in stripping mode (see Sec. 1 for a discussion of stripping mode; removing surfactant from Plateau borders such that they become leaner is an inherent part of stripping mode). Even when reflux is considered though, it is necessary to establish what is happening on the surface first, and then use those surface conditions to determine how surfactant is transported in the bulk.
On the surface though, as was the case with Vitasari et al. (2013b), we have already mentioned that there is a competition between Marangoni flow and film drainage. Due to the dominance, at least early on, of Marangoni flow, there is a flow on the surface in the direction from the Plateau border towards the centre of the film.
The physics that the model describes is as follows. In general a gradient in surfactant concentration is present on the surface. This gradient then is what sets up a flow field, and that flow field is what causes a convective flow also in the bulk. Comparatively close to the surface, the bulk flow carries fluid from the Plateau border towards the film. This fluid is richer in surfactant than the film itself, as reflux in foam fractionation tends to keep the Plateau border's concentration (c Pb ) higher than the bulk film concentration (c 0 ). Here in fact we work in a dimensionless system in which c Pb becomes unity, and c 0 also turns out (as we discuss later) to be the same as dimensionless Γ 0 (with Γ 0 < 1 here). Thus a new region in the bulk with the Plateau border's concentration appears, and is carried towards the centre of the film. However to compensate the incoming flow from the Plateau border, a leaner concentration fluid exits the film into the Plateau border. The latter flow is not from locations near the surface (z = δ), but instead mainly from locations closer to the midplane of the film (z = 0) with concentration c 0 as we have said.
As can be seen in Fig. 1b, the film is divided into two regions with distinct concentrations. One region has area A b . The other has area A t − A b , where A t here is the total area, or more specifically A t is the total area of a half-length and half-thickness of film. Owing to the way in which we make the system dimensionless, in fact A t is identical to δ as given by Eq. (12), or if the film is not draining, A t and δ are fixed at unity. In any case, once we know these areas A b and A t − A b we can also figure out how much surfactant is in the film. The boundary that divides these regions is to be called the separatrix, and to determine what the areas are we need to work out what the evolution of the separatrix is.

Calculating evolution of separatrix
As can be seen from Fig. 1a, the boundary which separates bulk and Plateau border concentrations is initially a vertical line which passes through x = 1 at the edge of the film (where it meets the Plateau border). As already alluded to, we call this boundary a separatrix. However during the foam fractionation process with reflux, the shape of this boundary changes continuously due to the effects of convective flow in the bulk. As mentioned, close to the surface, a flow of uniform concentration (in dimensionless form, c Pb ≡ 1) is pulled towards the centre of the film, while around the midplane of the film, fluid (with lower concentration c 0 ) moves out of the film towards the Plateau border. In the no film drainage case, the amount of fluid entering matches the amount leaving. However in the case with film drainage, the amount leaving is always slightly more than the amount entering. Nonetheless, since the fluid leaving tends to be leaner in surfactant than the fluid entering, there is still a possibility to use reflux to recover more surfactant in the film.
As mentioned previously, we need a model for how the separatrix evolves and then knowing the shape of the separatrix we must calculate the areas A b and A t − A b . As a result, we track a number of initially uniformly distributed material points on the separatrix with time. The general equations for how these material points evolve and hence how the separatrix is convected are as below (the details of implementing these equations numerically can be found in supplementary Sec. S 3): where x sep and z sep are x and z positions of the material points on the separatrix and u and w are velocity fields in the x and z directions. However at initial time, the velocity field Eq. (5) turns out to be singular, leading instantaneously to an infinite velocity (Bird et al., 1960). Hence, at early times, the numerics are difficult to handle and we need a bespoke method to evolve the separatrix early on. The relevant method is addressed in supplementary Sec. S 2.
Having defined the separatrix, it is now easy to calculate the size of the respective regions containing surfactant with the Plateau border's concentration and containing surfactant with the initial bulk concentration. This is discussed next.

Total amount of surfactant present in the foam film
In what follows we determine the size of the bulk regions containing surfactant with the Plateau border's concentration and containing surfactant with the initial bulk concentration, the regions themselves being sketched in Fig. 1b. The amount of surfactant in the bulk then immediately follows. However this does not represent the total amount of surfactant in the film, since we must also account for surfactant on the surface. This is again determined in what follows.
At a specific dimensionless time t, the total dimensionless amount of surfactant contained in the film, S T , is the sum of the surfactant S S on the film surface and S B in the film bulk In particular S S can be obtained from the following formula Here, S S is identical to what has previously been worked out in Vitasari et al. (2013b). However the dimensionless amount of surfactant S B contained in the film bulk is a new concept not considered by Vitasari et al. (2013b) and can be calculated as follows: where S is solubility parameter (Eq. (3)), c is dimensionless concentration of surfactant in the bulk. Note in particular the prefactor S appearing in this equation which follows owing to the way the system has been nondimensionalised. An analogous prefactor appears in the work of Rajabi and Grassia (2023).
Within Eq. (17), c has different values either side of the separatrix (in dimensionless variables, unity on one side, and c 0 on the other), but in each of those regions c itself is spatially uniform. Hence in each region, c can be taken outside the integration, and the integrals then merely compute areas, respectively A b and A t − A b . At the initial instant, A b = 0 and S B = S c 0 A t . Immediately after that, A b starts to grow and A t − A b falls. Note also the difference between c in this study and in Rajabi and Grassia (2023): c(x, t) in Rajabi and Grassia (2023) was in instantaneous equilibrium with Γ(x, t), but here it is not, remaining fixed instead at either unity or c 0 .
Thus far in this section we have considered surface S S and bulk S B contributions separately. However it is important also to understand how they are coupled. This is discussed in Secs. 2.5.1 and 2.5.2.

Extent of coupling between surface and bulk
Above we mentioned a difference between the present model and the work of Rajabi and Grassia (2023). Note another important difference from Rajabi and Grassia (2023) here. The flow fields that convect surfactant in the bulk here can be determined entirely from knowledge of Γ and δ (see Eqs. (5) and (6)). However Γ and δ both evolve in the present model entirely independently of solubility S as Eqs. (10) and (12) make clear. Consequently there is only one-way coupling here: the surface drives the bulk, but the bulk does actually not influence the surface. The separatrix shape that we compute is S independent, as is the integral term in Eq. (17). The only S dependence in S B therefore is due to the multiplicative prefactor S outside the integral. This means that we can solve for the separatrix shape just once, and we then know S B for all S values: in effect we are post-processing the results of Vitasari et al. (2013b) here. Of course the value of S S here is also independent of S and so must be the same S S as computed by Vitasari et al. (2013b) for an insoluble case. Since S T is nothing more than the sum of S S and S B , it can also be obtained for all S by doing a computation just once.
The situation is rather different from the regime considered by Rajabi and Grassia (2023). Coupling was much stronger there such that the bulk could also affect the surface. The evolution equation for Γ that resulted then depended explicitly on S, so to evaluate the various amounts of surfactant present (S S , S B and S T ) it was necessary to compute a separate solution for each different value for S, a rather more laborious process.

Relating bulk surfactant with surface surfactant
Returning to the problem at hand, thus far we have explained how to obtain the amount of surfactant in the bulk and the amount of surfactant on the surface, but not specifically how it might be possible to relate the two. In order to relate the surface and bulk surfactant concentrations, generally speaking we need an adsorption isotherm (Chang and Franses, 1995;Hiemenz and Rajagopalan, 1997;Shchukin et al., 2001). Especially when we are dealing with proteins rather than simple surfactant molecules, isotherms can take rather complicated forms (Fainerman et al., 2003;Fainerman and Miller, 2005;Gochev et al., 2021). Specifically what the isotherm does (Butler, 1932) is to relate the equilibrium amount of surfactant on the surface to the equilibrium amount of surfactant in the bulk.
However, in the dimensionless system with which we are working, the equilibrium amount of surfactant on the surface of the Plateau border and in the bulk of the Plateau border are both unity by construction. Any isotherm that we select must respect that. Even with that constraint though, there are still different isotherms that could be used. However following Rajabi and Grassia (2023), we will simplify the model and use what we call a global Henry adsorption isotherm. The global Henry isotherm in dimensionless form then requires that at equilibrium Γ = c. However a feature of the surfactant mass transfer model employed in the present work is that there is no general requirement at any instant for there to be equilibrium between surface and bulk. Equilibrium might still apply between the surface and a subsurface immediately adjacent to it. However diffusion is considered too slow on times scales of interest for equilibrium across the entirety of the bulk to be achieved.
14 Hence, in the specific model used here, we only ever utilise the isotherm to relate the amount of surfactant on the film surface initially with the amount in the film bulk initially. Using the global Henry isotherm this turns out to be in dimensionless form Γ 0 = c 0 , where Γ 0 and hence c 0 are necessarily less than unity: owing to reflux through the Plateau borders, the film starts off leaner in surfactant than the Plateau border. Any other isotherm could be chosen and would just give us a rather more complicated relation (Fainerman et al., 2003;Fainerman and Miller, 2005;Gochev et al., 2021) between Γ 0 and c 0 (see Sec. S 5 for an example). The requirement to have Γ 0 < 1 and c 0 < 1 in the dimensionless system here would however be retained.
By coupling the isotherm with the evolution with time of the separatrix shape, we now have a definitive formula for the amount of surfactant in the bulk, namely where the global Henry isotherm has been assumed, i.e. c 0 is the same as Γ 0 . We then use Eq. (16) to obtain S S , and Eq. (15) to obtain S T .
The way we proceed here is to set various different values for Γ 0 . The value of Γ at any given position and time depends of course on Γ 0 , but is, as we have mentioned, independent of the bulk, i.e. independent of S. Since flow fields depend on Γ (Eqs. (5) and (6)), and since flow fields also advect the separatrix (Eqs. (13) and (14) 15)) is affected by different Γ 0 is therefore less straightforward than determining how the evolution of S T is affected by different S.
Having now quantified the amount of surfactant, in the next section we present how we define recovery and enrichment of a fractionation process with respect to the parameters discussed earlier.

Recovery and enrichment
Recovery and enrichment are two quantities which are often used to evaluate the performance of a foam fractionation process (Blesken et al., 2020;Stevenson and Jameson, 2007). In this study, total amount of surfactant, S T is a measure at any instant of the recovery per foam film (or in fact part of a foam film, as we are using film half-length and half-thickness in our model). It also can be converted to the conventional recovery definition by specifying the number of foam films leaving the fractionation column (Rajabi and Grassia, 2023). Meanwhile, enrichment is the ratio between surfactant concentration in the foamate to the initial feed solution concentration and hence, enrichment can be quantified as S T /(S c 0 A t ), where S is solubility parameter (Eq. (3)), c 0 is initial solution concentration (and in our case is equal to Γ 0 ) and A t is instantaneous cross sectional area of a film half-length and half-thickness. We can see from the form of the equations (Eqs. (15), (16) and (18)), enrichment and recovery are dependent on solubility parameter, initial surfactant concentration and film half-thickness (given via Eq. (12)). Results for recovery and enrichment will be discussed later (Sec. 3.6).
Now, having defined the model, we solve it numerically using the procedure that we have already established in previous work (Rajabi and Grassia, 2023) along with some additional methodology to evolve the separatrix (see also supplementary Secs. S 2 and S 3) and the parameter values that we use are given in supplementary Tables S 2 and S 3 within Sec. S 4.1. Benchmarking is also done within supplementary Sec. S 4.2, so we turn now to results.

Results
In this section, results are discussed in the following order. We start by considering the evolution with time of the total amount of surfactant on the surface of a foam film in Sec. 3.1. Then Sec. 3.2 discusses the evolution within the film bulk of the so called separatrix with time, while Sec. 3.3 explains the evolution with time of A b /A t , where A b is the area of the region containing material initially in the Plateau border's bulk but then advected into the film. This is then normalised by the total cross sectional area A t . The effects of initial surface concentration Γ 0 and solubility parameter S on the evolution of total amount of surfactant in the film S T measured relative to initial amount of surfactant S T,0 are discussed in Sec. 3.4. We then analyse in Sec. 3.5 the effect of time evolution upon systems with various different initial surfactant concentrations and different solubility parameters. Then Sec. 3.6 considers recovery and enrichment in a foam fractionation process.

Evolution with time of the total amount of surfactant on the surface of a film
The evolution with time of the dimensionless amount of surfactant on the surface S S (Eq. (16)) is simulated and plotted in Fig. 2. More specifically we plot 1 − S S , this being a quantity which we know decays over time. The evolution of S S is the same as what happens for an insoluble surfactant and has previously been investigated by Vitasari et al. (2013b). However, due to the effect that evolution of surfactant on the surface has on the advection of surfactants in the bulk, we discuss S S here in brief. The results for the cases with no film drainage and with film drainage are displayed in Fig. 2a and Fig. 2b, respectively.
By comparing Fig. 2a with Fig. 2b, at early time, the evolution of S S is very similar for the cases not including or including film drainage effects. This is to be expected from Vitasari et al. (2013b), due to the dominance of Marangoni flow on the surface of the film over any film drainage effects at early times. However, as time continues to evolve and Marangoni flow becomes weaker, the competition between Marangoni flow and film drainage in the case with film drainage (Fig. 2b) slows down the evolution of S S slightly, compared to the no film drainage case (Fig. 2a). The main differences between the two cases are observed at later times though. In the case without film drainage (Fig. 2a), the surface evolves quickly towards a uniform surfactant distribution, in  Table S 3. dimensionless form Γ(x) → Γ Pb ≡ 1, with no surfactant concentration gradients remaining to keep pulling material onto the surface (hence S S = 1). On the other hand, when film drainage is considered, S S will only approach this final state rather more slowly. This is due to the fact that, at later times, a quasi-steady state between weak remaining Marangoni effects and slow film drainage is reached (Vitasari et al., 2013b). This prevents the surface from reaching a completely uniform concentration, at least as long as the film keeps draining. Fig. 2 also shows that the parameter Γ 0 affects at least slightly the dimensionless time it takes for the film surface to reach a uniform surfactant distribution without film drainage. Here, as Γ 0 is increased, leading to less discrepancy between the Plateau border surface and the film surface, the time to reach a uniform surfactant distribution on the surface is less, albeit this time is only a weak function of Γ 0 . On the other hand, in the case with film drainage it is apparent that the quasi-steady state once it is attained is independent of Γ 0 , as the curves for all the different Γ 0 values collapse together.

Evolution with time of the separatrix shapes
Since in this study surfactants are considered to be soluble, they are present in the bulk too. Thus, it is necessary to understand how surfactants in the bulk are transported due to the advection flow driven by surfactant transport on the surface. The results presented in Fig. 3 reflect the passive advection of points in the separatrix, where recall the separatrix is the boundary separating material that has arrived from the Plateau border from material that was originally in the film. These material points are initially distributed along the line at the edge of the film (x = 1), but move due to the advective flow in the bulk, generated either by Marangoni stresses on the surface alone, e.g. the no film drainage case (Fig. 3a for Γ 0 = 0.1), or due to the interplay between those Marangoni stresses and film drainage, e.g. the film drainage case (Fig. 3b again for Γ 0 = 0.1). As mentioned already, in this study, transport in the bulk is considered to be purely advective, because surfactant diffusive transport from the surface to the film bulk has been neglected (large Pe ∆ limit). Overall, it can be seen in Fig. 3 that the separatrix is pulled towards the the left, i.e. towards the centre of the film over time, at least at locations close to the surface (z = δ). However for locations close to the midplane of the film (z = 0), it is pulled to the right over time and ultimately out of the film. This is caused mainly due to the effect of Marangoni flow which is leftward on the surface, and rightward near the midplane. Film drainage if present, also competes with Marangoni near the surface, but cooperates with Marangoni near the midplane.
Results for the case with no film drainage, in Fig. 3a, show that, at early times, Marangoni flow is rapid, so the separatrix evolves quickly initially. At later times, Marangoni flow becomes much slower as the surfactant surface concentration gradients weaken as a result of the gradual enrichment of surfactant on the film surface. The separatrix shape therefore evolves increasingly slowly over time, as Fig. 3a shows, until it reaches a final steady shape for which Marangoni effects are no longer present because a surface uniformly covered in surfactant (Γ → Γ Pb ≡ 1) is achieved.
After that, there is no more convective exchange of material between the film and Plateau border.
In the case in which film drainage is included in the model along with Marangoni effects (Fig. 3b), at early times, flow due to Marangoni effects dominates film drainage. However as a result of an ongoing competition between Marangoni-driven flow and film drainage opposing it, the evolution of the separatrix shape is slightly slower.
Proceeding towards later times, the decrease in the vertical coordinate of the leftmost and topmost point in the separatrix (point at the film surface; see Fig. 3b) now shows the film becoming progressively thinner as time evolves, due to film drainage effects. This same leftmost and topmost separatrix point also of course migrates horizontally. However from about t ≈ 5 onward, the material on the surface of the film is barely moving horizontally at all, due to film drainage and Marangoni effects coming into a quasi-steady balance on the surface (Vitasari et al., 2013b).
Although the top left of the separatrix then no longer moves leftwards, material points lower down in the separatrix (within the film bulk) are still moving rightwards towards the Plateau border, as a result of ongoing film drainage in the bulk. As well as moving rightwards though, these same points are also moving downwards and this is the main effect we see in Fig. 3b at later times.
In Figs. 3c and 3d we show analogous data but for Γ 0 = 0.5 (instead of Γ 0 = 0.1). The main effect we see is that the separatrix is pulled less strongly to the left as Γ 0 increases. There is also a weak effect in the z direction (evident by comparing the right hand end of the separatrix in Figs. 3c and 3d with the right hand end of the separatrix in Figs. 3a and 3b). At any given time, increasing Γ 0 seems to move the right hand end of the separatrix downward very slightly relative to cases with smaller Γ 0 . However this is a much weaker effect than what is seen in the horizontal.
In summary, at early times, there is rapid surfactant exchange between the Plateau border and film bulks, dominated by Marangoni effects pulling material with dimensionless concentration c Pb ≡ 1 from the Plateau border into the film, at least for locations near the surface. Meanwhile for locations closer to the midplane of the film, material of dimensionless concentration c 0 (with c 0 = Γ 0 here) is pulled out of the film into the bulk of the Plateau border due to continuity. At later times, material inside the film is being pulled into the bulk of the Plateau border throughout, although locations near the midplane tend to be moving faster than those near the surface. As has been noted, different regions within the film have different concentrations, although those concentrations do not themselves evolve with time. As a result, by calculating just the areas of those regions we are able to calculate total amount of surfactant in the bulk of the film. Therefore in the next section we focus on these areas and how they evolve over time.

Evolution with time of area ratio A b /A t
In this section, the evolution with time of the area ratio A b /A t is presented. Recall that A b here is the area within the film containing material that was initially in the bulk of the Plateau border. Meanwhile (see Sec. 2.3), the total area A t has a straightforward evolution which is due to the film drainage: indeed in a no film drainage case, it is constant, equal to unity in the dimensionless system used here. Data are presented in Fig. 4. Specifically data for cases without film drainage and cases with film drainage are found in Fig. 4a and Fig. 4b, respectively. Results with no drainage, with A t now equal to unity at all times, are shown in Fig. 4a. We see that A b /A t increases significantly at early times. At later times though, A b /A t eventually reaches a final steady value. The final steady value of A b /A t is dependent on the dimensionless initial uniform surfactant concentration along the film surface, Γ 0 . Smaller values of Γ 0 (less surfactant initially on the film surface) will lead to higher values of A b /A t being achieved, i.e. more material being pulled from the Plateau border surface onto the film surface.
In the case where film drainage is included, A b /A t has a similar behaviour at early times (as observed in Fig. 4b), due to the dominance of Marangoni effects over film drainage effects early on. Indeed Marangoni effects manage to produce a significant amount of mass transfer even before the film has had time to drain substantially. However, as has also been observed in Secs. 3.1 and 3.2, the early time evolution is now slightly slower due to the competition between Marangoni flow and film drainage. Major differences in A b /A t for the two cases are however observed at late times, when the film drainage effect can no longer be neglected (see Fig. 4b).
To analyse the effect of Γ 0 on A b /A t , it is helpful now to write the rate of change of A b /A t with time as follows: As long as Marangoni flow dominates, d(A b /A t )/dt is always positive, so that A b /A t increases with time. At this stage, both terms on the right hand side of Eq. (19), i.e. both dA b /dt and −(A b /A t ) dA t /dt are positive. However, at late time for the case with film drainage, the second term on the right hand side is positive, but the first term on the right hand side can become negative. Hence, the sign of Eq. (19) depends on the relative magnitude of these two terms. Since in the cases with smaller Γ 0 , the value of A b /A t at any given time is larger (i.e. more surfactant exchange has taken place), the second term on the right hand side of Eq. (19) is larger and hence A b /A t remains 20 an increasing function of time. The opposite happens for the cases with larger Γ 0 , the value of A b /A t is then lower, and hence the second term on the right hand side of Eq. (19) is less important than the first term: the right hand side of Eq. (19) is then negative, and so the value of A b /A t can decrease at long times. Knowing the value of A b /A t and how it behaves with time is of interest, since A b /A t turns out to determine the average concentration of dissolved surfactant in the film bulk which can be obtained from c 0 + (c Pb − c 0 )A b /A t . Of course in the dimensionless system here, c Pb ≡ 1 and c 0 = Γ 0 .
In summary in the case with no film drainage, at early times, A b /A t increases quickly at first, increasing the overall amount of surfactant in the film, since material of concentration c 0 is being substituted by material of higher concentration, c Pb ≡ 1. At later times, the system reaches a steady state, where A b /A t reaches a final steady value. Furthermore, as we decrease Γ 0 , the earlytime rate of growth of A b /A t becomes faster, leading also to higher final values of A b /A t in the late-time limit. In the case where film drainage effects are present, a slightly slower initial increase is observed for A b /A t compared to the no film drainage case. At later times, once the film drainage effects start to become more significant in relative terms, two different situations can be observed depending on the value of Γ 0 . In cases with smaller Γ 0 , the value of A b /A t continues to increase, while it decreases at longer times in cases with larger Γ 0 .
Having now considered the areas of two regions within the bulk with uniform but different concentrations, now we can calculate the amount of surfactant in the bulk, as well as the overall amount of surfactant on the surface plus the bulk. This is discussed in the next section.

Total amount of surfactant relative to initial amount
Now, we proceed to analyse the influence of the solubility parameter S (Eq. (3)) and initial surfactant concentration Γ 0 , on the total amount of surfactant contained in the film S T (surface plus bulk). Specifically, we plot total amount of surfactant (S T ) relative to its initial value (S T,0 ), for different S and Γ 0 (see Figs. 5 and 6). This ratio S T /S T,0 is a measure of how reflux affects the recovery of a fractionation process over time relative to the initial recovery prior to any benefit from reflux. The larger this ratio, the more the system benefits from reflux. However the other relevant factor in cases with film drainage is thinning of the film, which can result in a decrease of the total amount of surfactant in the film, and subsequently, a decrease in the recovery over time. First, we present the results of the no film drainage case in Sec. 3.4.1. Then, the case with film drainage will be discussed in Sec. 3.4.2.
In this study, using relevant data from literature (see Tables S 2 and S 3), a base case estimate for the solubility parameter S (as per the definition in Eq. (3)) turned out to be 3.09. To investigate the effect of solubility parameter being varied about the base case, we selected values S = 3 and S = 30 for the no film drainage case. In this particular case, we did not explore solubility parameters smaller than S = 3, because the completely insoluble case has been solved by Vitasari et al. (2013b). We know already what happens there: S T /S T,0 grows over time from unity to Γ −1 0 , 21 at least in the absence of film drainage. In the case with film drainage, behaviour is more complex, and so to elucidate the behaviour a little more, we look at a wider domain of S, namely S = 0.03, S = 0.3, S = 3 and S = 30.

Case without film drainage
The results of the no film drainage case are presented in Fig. 5 for two different solubility parameters, 3 and 30 as we have said. What we see in each case is that S T /S T,0 increases with time and then eventually reaches a final steady state value. The final steady state value becomes larger as the initial surfactant amount Γ 0 becomes smaller. Moreover the final steady state value also becomes larger as the solubility S decreases. However the final steady state value of S T /S T,0 always falls short of Γ −1 0 (the value attained for an insoluble surfactant (Vitasari et al., 2013b)). The final amount of surfactant in the film can be calculated by considering the final values resulting from Eqs. (15) to (18). This then leads to the equation where A b now specifically denotes a final area. This final S T value can also be written in the form . In our dimensionless system Γ Pb and c Pb are unity, A t is also unity (in the absence of film drainage), and c 0 is the same as Γ 0 . Hence, with those substitutions Recall also (see Sec. 2.5.1) that the only S dependence here is the dependence showing explicitly: the value of A b is sensitive to Γ 0 but not sensitive to S.
By similar arguments, in our dimensionless system S T,0 = (1 + S)Γ 0 . It is clear now that the ratio S T /S T,0 always falls short of Γ −1 0 , and the amount it falls short grows as S grows. In a hypothetical case in which A b /A t became as large as unity, S T /S T,0 could attain the value Γ −1 0 . However we know from Fig. 4a that A b /A t never becomes that large. The other way to prevent S T /S T,0 falling short of Γ −1 0 is to have Γ 0 itself approaching unity. However the difference between the initial surfactant concentration on the film and the concentration in the Plateau border is then so small, that there is essentially no benefit to derive from reflux.
In summary, the total amount of surfactant in the film S T increases when solubility S, is increased, but the relative amount, S T /S T 0 actually decreases. We observe that high solubility results in having more surfactant in the film initially, but simultaneously that reduces the factor by which reflux can then increase the recovery over time.

Case with film drainage
Regarding the case in which film drainage is included (Fig. 6), at early times the behaviour of S T /S T,0 is similar to what is seen when film drainage effects are neglected. However the rate of increase of S T /S T,0 is slightly less than in the no film drainage case, because of the competition between Marangoni flow and film drainage, which slightly slows down the evolution of S T .
At later times though, a decrease in S T /S T,0 is observed. This is due to the fact that, after reaching 22 Transport of convected soluble surfactants on and within the foam film in the context of a foam fractionation process  Table S 3. a (quasi-steady) balance on the surface, the main effect influencing the amount of surfactant is film drainage from the bulk. Hence, the total amount of surfactant or equivalently surfactant recovery of the film decreases from this time onward. For the less soluble surfactant cases (Figs. 6a and 6b), the late time decrease in S T /S T,0 is comparatively small. This is due to the fact that in cases with only moderate solubility, although film drainage is indeed removing material from the bulk, there is now barely any surfactant in the bulk to start with. However the reduction at late times is more significant in the more soluble surfactant cases. The other important factor affecting the late-time behaviour is the initial surfactant concentration. In fact, most of the surfactant loss from the bulk is from locations near the midplane of the film containing surfactant with the initial bulk concentration c 0 (equal to Γ 0 in the dimensionless system here). As a result, films with lower initial bulk concentration tend to lose only surfactant lean material, whereas films with higher initial concentration stand to lose more surfactant.
Once sufficient material is lost from the bulk, the surfactant recovery S T for a film might even become less than the initial amount S T,0 . As Fig. 6 shows, this typically happens for the cases with higher solubilities S and higher initial surface concentrations Γ 0 . Even when that happens though, there might still be a benefit for the film in enrichment terms, as we will see later on.
In summary, in cases with film drainage, surfactant recovery reaches a maximum around the time that a quasi-steady state is reached, then decreases. In cases with higher solubility parameters, due to the same reason mentioned in Sec. 3.4.1, the increase of S T /S T,0 up to the quasi-steady state is less. However its decrease from that time onward is higher: when S is large, most of the surfactant is contained in the bulk, and we are specifically losing surfactant from the bulk. Overall, more soluble cases are more sensitive to film drainage.

Effect of time evolution in cases with film drainage
In this section, we investigate in more detail how S T /S T,0 and also S T itself are affected by the fractionation process time evolution. We focus on the case with film drainage, since this is the case which exhibits interesting non-monotonic behaviour. We have chosen just two values of the solubility parameter S, namely 3 and 30, as these are cases in which the non-monotonic behaviour is more evident. The times when S T /S T,0 and S T reach a maximum for various values of Γ 0 are identified, along with the maximum values themselves. Moving to later times, values of S T /S T,0 at 10 and 100 time units are also presented for selected S and Γ 0 within Tables 1 and 2.
Looking at Tables 1 and 2 separately, reveals the fact that in each case there is a certain time corresponding to the maximum S T or equivalently maximum S T /S T,0 , and this time is highly dependent on the value of Γ 0 . For lower Γ 0 values, longer time is required for the maximum to happen. What leads to the maximum in S T /S T,0 is the following. Initially, Marangoni effects dominate, hence S T /S T,0 increases. However film drainage is also taking surfactant out from the bulk. After some time, the amount of surfactant being removed from the film due to film drainage  Moreover, it is confirmed that, as Fig. 6 also shows, systems with lower Γ 0 eventually reach a higher maximum value of S T /S T,0 . This means that the cases with lower Γ 0 gain in relative terms more surfactant due to reflux, despite the fact that total amount of surfactant in these cases is less than cases with higher initial concentrations.
Looking now towards even longer times, e.g. t = 10 and t = 100, less of the surfactant that was gained is subsequently lost in the cases with lower Γ 0 . This is again due to having leaner bulk in these cases. For instance, in the case with Γ 0 = 0.1 and solubility parameter S equal to 3, at 10 and 100 times units respectively, S T /S T,0 is approximately 1% and 15% less, relative to its maximum amount. However the analogous amount for Γ 0 = 0.9 is around 8% and 36% less at 10 and 100 time units, respectively. Leaner solutions therefore not only benefit more from reflux up to the maximum, but also lose less surfactant as the process continues. As a result enrichment (see definition in Sec. 2.6) in these cases can be high without much decrease over time in their recovery (to be discussed later on in Sec. 3.6). On the other hand, solutions that were richer initially e.g. Γ 0 = 0.9 lose more surfactant later on, and indeed they can lose so much surfactant that they end 25 up by t = 10 or t = 100 with less than they had initially. Overall, it can be seen that leaner solutions benefit more from foam fractionation with reflux. Now, instead of looking at Table 1 and Table 2 individually, we can compare them to gauge the effect of different solubilities. In fact the initial total surfactant amount S T,0 , in the systems which contain higher solubility surfactants are significantly higher. Hence, as far as recovery is concerned, even without benefitting from reflux, surfactants with higher solubilities can be recovered to a greater extent than lower solubility surfactants. Reflux merely increases that recovery, i.e. the maximum S T is even higher than S T,0 . However, as can be seen in Table 2, the time when S T /S T,0 for more soluble surfactants reaches a maximum is significantly less than cases with lower solubilities. This is due to the fact that in these cases the contribution from the film bulk is more important than in low solubility cases. As a result, the effect of losing surfactant from the bulk due to the film drainage can be seen sooner.
This is more obvious in systems with higher Γ 0 which need very little time at all to reach a maximum in S T /S T,0 . For instance, in the case with Γ 0 = 0.9 and solubility parameter S equal to 30 say, the time required to reach the maximum is just 0.01 time units and the maximum amount of S T itself is only around 0.02% higher than the amount of surfactant at initial time S T,0 . In addition, in higher solubility cases, more surfactant is lost at longer times. For instance, approximately 12% and 49% of surfactant relative to its peak amount will be lost after 10 and 100 time units in this last mentioned case (Γ 0 = 0.9 and S = 30). Thus, for these cases there is barely any benefit for having long residence time if increasing the recovery is our primary goal in foam fractionation. However, as is discussed next (Sec. 3.6), reduction of recovery even in these cases can still be offset by increased enrichment as the film thins.

Recovery and enrichment
We now consider recovery and enrichment which as was discussed in Sec. 2.6 are important parameters for assessing performance of a fractionation process: formal definitions are given in Sec. 2.6. In Fig. 7, we have plotted recovery versus enrichment for various solubility parameters (0.1 ≤ S ≤ 50) at different times, assuming film drainage is present, and for two different initial surfactant concentrations (Γ 0 ), namely 0.1 and 0.5. We do not consider Γ 0 = 0.9 here, because reflux is expected to have only very modest benefit in that case.
We can see that there is a trade off between recovery and enrichment. As a result, to increase one, we generally need to decrease the other. High solubility parameter S leads to high recovery but low enrichment, whereas low solubility parameter leads to low recovery but high enrichment. Nonetheless, in the presence of reflux, recovery and enrichment benefit as time evolves: the recovery-enrichment curves are pushed upwards and to the right. However the extent of that benefit depends on the initial surface concentration (Γ 0 ): there is much more benefit when Γ 0 is small (Fig. 7a) than when it is larger (Fig. 7b).
The other observation is that in all cases, after around t = 10 further benefits in the recovery versus enrichment plot are limited: the t = 10 and t = 100 curves lie close to one another, so are a little difficult to distinguish in the plots. Between t = 10 and t = 100, as film drainage continues to remove material from the bulk, in the high solubility part of the plot there is a modest decrease in recovery, while in the low solubility part of the plot there is a modest increase in enrichment. To summarise then, the effect seen at longer times, is an enhancement of enrichment at the expense of recovery. Long times (i.e. many units of dimensionless time) are needed to see this though, since (as already noted in Sec. 2.1) drainage flow is comparatively slow.
It is worth remembering also here that throughout data have expressed in terms of dimensionless time. The conversion back to dimensional time using typical parameter estimates is given in Table S 2. For the given parameter values, dimensionless t = 100 corresponds to only around 1 s of dimensional time. Even this is short compared with a typical residence time that might be expected for a foam film in a fractionation column. Hence not just in Fig. 7, but in Figs. 2 to 6 and in Tables 1  and 2 also, it is actually the long-time limiting behaviour that is of particular interest for typical fractionation applications.

Conclusion
Simulation of convected soluble surfactant transport on and within a foam film has been carried out in the context of a foam fractionation process with reflux (or analogously, a foam fractionation process operated in stripping mode). This study builds on related models developed and published earlier (Rajabi and Grassia, 2023;Vitasari et al., 2013b). However the main novelty of the present study is that here mass transfer of soluble surfactant in the bulk of a foam film due to convective flow is taken into account. The convective flow happens as a result of an interplay between Marangoni flow and film drainage on the surface of the foam film, albeit with the resulting flow fields then extending into the bulk. The parameters used for the simulation purposes were taken from relevant literature (Fainerman et al., 2020;Gochev et al., 2013Gochev et al., , 2021Lexis and Willenbacher, 2014;Miller et al., 2004;Pradines et al., 2009): see supplementary Sec. S 4.1 for details. Here also dimensionless parameters that were identified as being particularly relevant to system behaviour were the solubility parameter S and initial surface surfactant concentration Γ 0 , so various values of these parameters have been considered. Note that we are referring to "surfactant" generically in this work, but the main application considered here is protein separation, e.g. β-LG.
In general, this study tackles a specific limit, namely large Pe ∆, in which diffusion of surfactant from the surface of a foam film into the bulk is ignored, so convection dominates. For simplicity and following on from the work Rajabi and Grassia (2023), we have selected a Henry isotherm to relate the bulk and surface concentrations. However we could have used another adsorption isotherm at the expense of making the relationship between bulk and surface non-linear (Fainerman and Lucassen-Reynders, 2002;Fainerman et al., 2003Fainerman et al., , 1996: adsorption isotherms for proteins in particular are often rather complex (Fainerman et al., 2003;Fainerman and Miller, 2005;Gochev et al., 2021). Here though specifically the isotherm is needed to relate initial surface and initial bulk concentrations. The consequence of switching to a non-linear isotherm would be to have less bulk concentration for a given surface concentration (see Sec. S 5 for an example). Thus we would recover less surfactant from the bulk than the model here predicts.
In addition, in line with previous literature (Rajabi and Grassia, 2023;Vitasari et al., 2013b), two different cases during this study have been considered, namely, a no film drainage case and a case with film drainage. The main difference between these two cases is that there is more surfactant on the surface and in the bulk of the film in the case without film drainage than with film drainage. Indeed, in the film drainage case, so much surfactant can be lost at late times, that even less is recovered than was present initially. Even so, there can still be a benefit in terms of enrichment.
In this study, within the bulk, we track a locus that we call a separatrix, which forms a boundary between two regions with uniform, but different concentrations. These regions are the one with the Plateau border's concentration and the one with the initial film concentration (the Plateau border's concentration is higher than that of the film, due to reflux). The separatrix evolves with time as a result of the convective flow extending into the bulk, thereby defines how the areas covered by the aforementioned regions change with time. Knowing the areas enables us in turn to calculate the amount of surfactant in the bulk and subsequently, the overall amount of surfactant in the film (surface and bulk together).
As surfactant evolution on the surface has previously been investigated elsewhere (Vitasari et al., 2013b), we mostly focus on what happens in the bulk and its effect on the overall amount of surfactant. It has been found that over time, the top of the separatrix is pulled towards the centre of the film, while for locations closer to the midplane of the film, the separatrix is pulled instead towards the Plateau border. This is as expected consistent with the direction of convective flow in the bulk. The parameter that has most influence upon the evolution of the separatrix is found to be the initial surfactant surface concentration Γ 0 that causes the driving force for setting up Marangoni flows. Lower surfactant surface concentrations create a larger driving force on the surface, hence, a stronger convective flow towards the centre of the film, at least for locations close to the surface. This then causes the separatrix to displace more and results in having more fluid with the Plateau border's concentration arriving in the bulk of the film.
It has also been found that higher initial surfactant concentration Γ 0 and higher solubility parameter S lead to having more initial overall amount of surfactant on and within the film, but less overall amount of surfactant relative to the initial amount. In the no film drainage case, the overall amount of surfactant relative to the initial amount increases up to a point and then reaches a steady state. However, in the case with film drainage, surfactant amount relative to initial amount decreases after reaching a maximum.
To analyse the effects due to drainage at later time when this decrease happens, we determined the time at which the overall amount of surfactant relative to the initial amount reaches a maximum, and also the decrease in amount of surfactant afterwards. It has been observed that in the cases with higher initial surfactant concentration Γ 0 and higher solubility parameter S, the time to reach maximum is less, and the amount of surfactant being lost from the film afterwards is greater. Thus, in cases with higher Γ 0 and higher S, the residence time optimising the surfactant recovery is much lower, and the negative effects of extending residence time to longer times are also more significant.
Hence we can say that these cases are more sensitive to the residence time (or equivalently to the fractionation column length). However the system has been made dimensionless on a scale such that typical residence times will likely be well beyond that optimum time. Moreover this is merely an optimum for recovery: longer residence times can still lead to better enrichment. Indeed, this study revealed that higher solubilities are more beneficial for surfactant recovery, and lower solubilities are more beneficial for enrichment of the foamate in a fractionation process. Systems with lower initial concentrations benefited more from reflux in terms of both recovery of surfactant and enrichment.
In summary, previous work (Vitasari et al., 2013b) looked at the insoluble case. In Rajabi and Grassia (2023), solubility was taken into account, but in the limit of small Pe ∆ (rapid surfactant diffusion across films, applicable for relatively small surfactant molecules and very thin films).
Here, we worked out another relevant limit, namely soluble surfactant with large Pe ∆ (convectiondominated limit, for bulky surfactant molecules within films that are not quite so thin). However in reality systems will presumably be in between these two limits. Moreover, just as Rajabi and Grassia (2023); Vitasari et al. (2013b) did, we simplified the film's shape, using instead the assumption of a flat film. As a result, even though we captured Marangoni and film drainage effects to an extent, the flow fields we determined are simplified ones. The adsorption isotherm has also been simplified as has been mentioned. Furthermore Plateau borders are treated as if they were surfactant reservoirs. However, despite the simplifying assumptions and the particular limit assumed, we were able to obtain predictions of fractionation performance, which can help to guide the separation by fractionation of surface-active materials, such as proteins.
There is an important point here however for a practitioner aiming to use fractionation to separate a big bulky molecule like a protein (likely to have a large Pe ∆ and hence involve convection-dominated transport): relying on convective transport to enrich a foam film significantly could well be challenging. Indeed convection is only able to enrich a bulk foam film to the extent that bulk liquid is transferred between a film and a surrounding surfactant rich Plateau border. Since the amount of liquid convected is always less than the total amount of liquid in the film's bulk, the extent of enrichment could therefore be rather limited. A less bulky molecule (e.g. a smaller surfactant molecule, instead of a protein) with a small Pe ∆ can be transported into the bulk of foam films via diffusion, and might therefore become enriched more readily.
The predictions however come with a number of caveats, that are important to mention. Although we have considered only convective surfactant transport in the film bulk, in reality Pe ∆ is large but finite (on the order of 1000 or so according to Table S 3). This means weak diffusion of surfactant will actually be present concurrently with the convection. That diffusion will of course be active where concentration gradients are largest. Clearly there are large concentration gradients at the separatrix itself, so a sharp concentration front there will be smeared out by diffusion even as the separatrix advects. Diffusion will also be present near the film surface. Material in the subsurface (immediately below the surface) can be in equilibrium with the surface, but is not in equilibrium with the remainder of the bulk slightly further from the surface.
However weak diffusion (whether near the separatrix or near the surface) can also couple with the convective flow fields. If diffusion causes surfactant to diffuse a small distance away from what is nominally a sharp front, then spatiotemporal differences in velocity will exist between the actual position of the surfactant and the nominal position of the front. Depending on the velocity gradients present in the flow field, this potentially could advect the surfactant even further away from the nominal position of the front.
Another caveat is that despite some similarities between the insoluble surfactant case of Vitasari et al. (2013b) and the large Pe ∆ case considered here (i.e. in both cases limitations are placed on surfactant exchange between surface and bulk), there are also differences. In the insoluble case, i.e. small values of S, even if surfactant did manage to transfer from surface to bulk, according to Eq. (17) there would never be much of it in the bulk, even if surface and bulk were to equilibrate. In the present situation, i.e. large values of Pe ∆, it is merely the case that Eq. (9) predicts surfactant diffuses off or onto the surface comparatively slowly, certainly on time scales longer than any Marangoni or film drainage flows might proceed. After several units of dimensionless time then, the separatrix will have reached a steady configuration (or in the case when film drainage is present, a quasi-steady configuration), and convection will have have come almost to a stop or at least will have slowed significantly. Even after that though, diffusion is still ongoing: provided S is not too small, significant amounts of surfactant can still eventually manage to enter the bulk diffusively, albeit as we have said, this occurs slowly compared to the prior convective transport. The surfactant distribution that results at the end of the convective transport stage is an initial condition for the diffusive transport stage.
In fact based on Eq. (9) we would need dimensionless times on the order of several multiples of Pe ∆ to see significant diffusion. A single unit of dimensionless time in the present problem, when converted back to dimensional variables (see Table S 2) is however only around 0.01 s. This is short compared to a typical residence time expected for a foam film rising through a fractionation column in a laboratory experiment, which could be on the order of seconds or even more. Hence the end of the convective stage is reached quickly, at least on a laboratory time scale. The dimensionless time scale for diffusion meanwhile is, as we have said, several multiples of Pe ∆. With Pe ∆ values on the order of 1000 as given by Table S 3 this corresponds to a time scale up to the order of a minute. On that sort of time scale, diffusive transport cannot be ignored, and the model already considered by Rajabi and Grassia (2023) (which incorporates significant diffusive surfactant transfer effects right from the outset) arguably then becomes more pertinent than the convective transport dominated model that has been considered here.

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.