A model for foam fractionation with spatially varying bubble size

A model is developed for a foam fractionation process, supposing that bubble size increases moving up the foam column. Predictions are obtained for liquid ﬂux through the foam column, liquid fraction at the top and enrichment of surface actives within the foamate. Compared to previous models with ﬁxed bubble size, fractionation performance is shown to be improved. A synergy is revealed between the eﬀect of the foam becoming drier (and hence richer in surface actives) as the foam column becomes taller, and the eﬀect of bubble size increasing moving upwards (which makes the foam even drier and hence richer still). A dimensionless parameter 𝑚 ∗ quantiﬁes the relative variation in bubble size between top and bottom for a foam column of typical height. Even small 𝑚 ∗ values lead to behaviour qualitatively diﬀerent from a system with no bubble size variation whatsoever, whilst increasing 𝑚 ∗ leads to even better fractionation performance.


Introduction
Foam fractionation is a promising separation technique (Lemlich, 1968(Lemlich, , 1972Stevenson, 2012Stevenson, , 2014 in many applications involving surface active materials. Such materials adsorb to bubble surfaces. Hence when a foam is created and the resulting foamate is collected (Martin et al., 2010), it is richer in surface actives than the original feed liquid was.
Improving the design and operation of fractionation systems can be facilitated by having physical models available (Du et al., 2000; Steven-* Corresponding author.
E-mail address: paul.grassia@strath.ac.uk (P. Grassia). son and Jameson, 2007;Hutzler et al., 2013;Tobin et al., 2014;Keshavarzi et al., 2022b). Unsurprisingly, given that fractionation involves adsorption of surface actives, models must necessarily incorporate a description of the physical chemistry of adsorption. In fact the target materials mentioned earlier for separation (proteins and biomolecules) can exhibit particularly complex adsorption behaviour because they can adsorb in many different states (Fainerman et al., 2003;Miller et al., 2004;Gochev et al., 2013). Over and above describing adsorption, additional physics is required within models. Injecting air into a liquid to form foam (as happens during fractionation) creates a multiphase system, and hence involves aspects of multiphase flow, which must therefore be incorporated within models also. A typical fractionation column contains a column of foam (hereafter "foam column") towards the top, and a bubbly liquid phase underneath: note the distinction made here between the "fractionation column" and the "foam column", with the latter filling just part of the former. Complex multiphase flow processes occur, as described below, both in the underlying bubbly liquid phase (Sarhan et al., 2017a(Sarhan et al., ,b, 2018b and in the foam itself (Verbist et al., 1996;Weaire et al., 1997;Koehler et al., 1999Koehler et al., , 2000. For instance, air can be introduced at different rates, and the resulting flow velocities within the bubbly liquid phase determine the balance between bubble break up and bubble coalescence events (Sarhan et al., 2017a(Sarhan et al., ,b, 2018b. This then determines the bubble size entering the https://doi.org/10. 1016/j.ces.2023.119163 Received 31 May 2023; Received in revised form 2 August 2023; Accepted 5 August 2023 foam and ultimately the surface area flux (Sarhan et al., 2018a) within it. It is in fact the foam rather than the bubbly liquid which is the focus of the present work, so the foam is discussed next.
Coalescence can also continue within the foam itself. However the extent of coalescence there can be sensitive to the amount of surface active material present (Keshavarzi et al., 2022b). On the film scale, this then impacts upon transport of surface actives along a foam film, which is itself a complicated process (Yeo et al., 2001;Vitasari et al., 2013;Rajabi and Grassia, 2023). However, in the context of fractionation, what is relevant is that typically having more surface active material present implies better foam stability and hence less coalescence (Keshavarzi et al., 2022b). Remember though that fractionation may often be used for quite dilute systems (King, 1980;Li et al., 2016;Tian et al., 2018) with comparatively little surface active material actually present.
Another parameter that is found to affect the fractionation performance (Stevenson, 2012(Stevenson, , 2014 is the foam liquid fraction. The drier the foam, the more enriched the foamate can be (Hutzler et al., 2013;Tobin et al., 2014;Grassia, 2023), because a specified amount of adsorbed surface active that is carried on bubble surfaces is then accompanied by very little bulk liquid.
To summarise then, knowing how liquid is distributed within and transported through the foam phase (Neethling et al., 2000) is therefore important for modelling fractionation performance. It turns out that such transport is governed by foam drainage theory, involving a balance between liquid convection by air, gravity-driven liquid drainage and capillary suction effects (Verbist et al., 1996;Weaire et al., 1997;Koehler et al., 1999Koehler et al., , 2000Grassia et al., 2001;Lorenceau et al., 2009). Foam drainage is typically an unsteady state process (Verbist et al., 1996;Cox et al., 2000;Neethling et al., 2005;Grassia et al., 2006;Brito-Parada et al., 2013), and a model for foam fractionation can be set up (Keshavarzi et al., 2022b) to capture that unsteady behaviour.
Nevertheless even simplified drainage theories (Grassia et al., 2001) (assuming e.g. liquid fluxes that are spatially uniform and also temporally steady (Hutzler et al., 2013;Tobin et al., 2014) or quasi-steady (Grassia, 2023)) can make useful predictions about fractionation performance: indeed quasi-steady state drainage behaviour tends to set up on a time scale that is short compared to the duration of a typical fractionation process (Grassia, 2023). It can then be shown for instance (Hutzler et al., 2013;Tobin et al., 2014;Grassia, 2023) that making the foam column taller also makes it drier and hence enriched in surface actives, whilst decreasing air flow velocity has the same effect.
The effect of increasing bubble size entering the foam is less straightforward. On the one hand, larger bubbles have lower specific surface area  and hence have less adsorbed surface active per total foam volume. On the other hand, larger bubbles drain more readily (Hilgenfeldt et al., 2001) and as a result have lower liquid content per total foam volume. In enrichment terms, it turns out that the latter effect outweighs the former (Grassia, 2023).
There was however one significant element missing from the simple models discussed above. Even though the effect of having various different bubble sizes was considered (Grassia, 2023), the bubble size was assumed to be spatially uniform throughout the foam column. In fact it has been found (Keshavarzi et al., 2022b;Tong et al., 2011) that bubble size actually grows moving up through a foam column, with the change in bubble size being more significant when the concentration of surface actives is low (Keshavarzi et al., 2022b). There is potential now for the spatially varying bubble size to impact fractionation performance, and thus likewise for foam stability to impact performance (Neethling and Brito-Parada, 2018).
The aim of the current study therefore is to incorporate spatially varying bubble size as considered by Keshavarzi et al. (2022b) into the comparatively simple foam drainage models considered by Hutzler et al. (2013); Tobin et al. (2014); Grassia (2023), and consequently to make novel predictions about fractionation performance. As we will see this leads, not just to quantitative differences in predicted performance, but qualitative differences as well. Moreover, as will be explained later, a synergy will be shown to arise between the effects of increasing foam column height and increasing bubble size.
The rest of this paper is organised as follows. Section 2 presents the models to be employed. Section 3 analyses the models, and then section 4 presents results. Finally section 5 offers conclusions. Additional information about the models is also provided in supplementary material.

Foam fractionation model
The simple model for foam fractionation used here involves three elements. The first element (section 2.1) is a description of liquid flux through the foam column, and the associated variation in liquid fraction over the height of the foam column. The second element (section 2.2) is the variation in bubble size. The third element (section 2.3) is the enrichment of the foamate, or equivalently the effective concentration of surface actives in the foamate. Here, as in Grassia (2023), the model will be presented in dimensionless form. Dimensional analogues are presented in sections S1-S2 in supplementary material. To give a sense of scale in what follows, we also indicate in due course what each dimensionless unit typically represents in actual physical units.

Model for liquid flux through and liquid fraction within the foam
We follow the same model for liquid flux through the foam column thru as was used by Grassia (2023), just with a minor modification to account for variation in bubble size. Models developed by Hutzler et al. (2013); Tobin et al. (2014) are similar, but the notation used in the model of Grassia (2023) lends itself more readily to considering variation in bubble size. The model can be written thru = air − 2 2 − −1 1∕2 ( 2 )∕ . (1) The terms on the right hand side represent respectively convection of liquid by the air, gravity drainage and the effect of capillary suction.
Here is foam liquid fraction. Also air is dimensionless air velocity (a typical value according to Grassia (2023) is 0.00195; see Table 1). Note that (see Grassia (2023) and also section S1 in supplementary material) one unit of dimensionless velocity corresponds to a physical velocity of around 0.056 m s −1 . As section S1 explains in detail, this velocity unit is essentially ∕ ( being surface tension and being liquid viscosity), with some additional prefactors associated with foam geometry. Section S1 explains that changing from one fractionation system to another has barely any impact on the value of this characteristic unit of velocity.
In addition is dimensionless bubble radius (discussed further in section 2.2) and is dimensionless vertical coordinate (measured upwards). The reason for using partial derivative notation with respect to is explained in section S1. Equation (1) applies over the domain 0 ≤ ≤ where denotes the dimensionless foam column height. Here one unit of dimensionless length (see Grassia (2023) and also section S1) corresponds to a physical length around 1.3 ×10 −3 m. Section S1 explains that this length unit is essentially ( ∕( )) 1∕2 (with being surface tension, being density of liquid in the foam and being gravity acceleration), but with some prefactors associated with foam geometry. Again section S1 explains that changing from one fractionation system to another has barely any impact on the value of this characteristic unit of length.
Note also in equation (1) the derivative that is acting on 2 rather than just on . Here 2 is in effect a dimensionless measure of the Plateau border cross-sectional area. Indeed the capillary suction pressure associated with Plateau borders (Weaire and Hutzler, 1999) is a function of the cross-sectional area rather than solely a function of liquid fraction. If is spatially uniform, then the model of Grassia (2023) is recovered. On the other hand, if varies spatially even at fixed , then the capillary pressure also varies, and equation (1) accounts for that. Another minor point we mention (Hutzler et al., 2013;Tobin et al., 2014;Grassia, 2023) is that equation (1) formally assumes a dry limit ≪ 1. What this means is that liquid volume fraction might vary in relative terms by orders or magnitude across the height of the foam column, whereas the relative variation of air volume fraction 1 − is less. Under these circumstances, air flux and air velocity are essentially the same, and thus air velocity air can be treated as spatially uniform. The forms of the gravity drainage and capillary suction terms within equation (1) also assume a dry limit. These terms would need to take more complicated forms in wetter foams Höhler et al., 2021).
Equation (1) needs to be solved with boundary conditions. At the bottom, we impose a condition = bot where we set bot = 0.36 (Grassia, 2023;Cantat et al., 2013). This boundary condition is used (Grassia, 2023) despite it extrapolating equation (1) beyond the dry limit. Typically however foam columns are tall enough so as to be dry over most of their height, away from the bottom.
At the top we impose a "no slip" condition that the velocity of liquid and air are the same. In other words thru = top air or equivalently top = thru ∕ air where top is the liquid fraction at the top. This supposes (see Grassia (2023) for more explanation) that (in line with experimental observations on a particular fractionation system by Keshavarzi et al. (2022b)) there is no significant bursting of foam films right at the top: this is what then requires liquid and air to have the same velocity there. Analogous boundary conditions also apply in the gravity thickening of suspensions (for detail, see Grassia (2023) and references therein, e.g. Fitch (1966); Usher and Scales (2005)). Returning to consider the case of foam fractionation, internal coalescence of bubbles (leading to spatial variation of bubble size) is still permitted (Keshavarzi et al., 2022b), and will be discussed in section 2.2.
What is apparent is that two boundary conditions have been imposed here, but equation (1) is only a first order differential equation. This implies that there must be a relationship between the foam column height and the liquid fraction at the top top . Establishing this relationship is an essential part of solving for the fractionation performance. The relationship in question is sensitive to spatial variation in the bubble radius, so this is considered next.

Model for changes in bubble radius
The novel aspect taken into account in the model employed here is spatial variation of bubble size resulting from bubble-bubble coalescence inside the foam. The work of Keshavarzi et al. (2022b) considered bubble size to vary linearly with vertical coordinate. We adopt the same approach here. Specifically we assume = 0 (1 + * ∕ base ). (2) Here 0 is the bubble radius entering the bottom of the foam column (a typical value according to Grassia (2023) is 0.25 dimensionless units; see Table 1). In addition base is a typical base case height of the foam column during a fractionation experiment (a suitable value according to Grassia (2023) is 40 dimensionless units; see Table 1): note that base is merely a typical value, and the actual foam column height might well differ from it. Recall also (see Grassia (2023) and section S1) that one unit of dimensionless length corresponds to a physical length around 1.3 × 10 −3 m. Finally * is a parameter that measures the stability of the bubbles against coalescence. The value of * depends on the amount of surface active agent present: higher concentration of surface actives implies smaller * , whereas lower concentration of surface actives implies larger * . Of course the concentration of surface actives must never be so low that the foam collapses entirely before it manages flow out of the fractionation system (Tobin et al., 2014): in order to enrich the surface active at all it is essential to collect at least some foamate.
Experimental data suggest that * might be in the domain 0.5 to 2 (Keshavarzi et al., 2022b), although simulations have considered values of * covering a domain 0.2 to 10 ( Keshavarzi et al., 2022b). Remember here that larger values of * would typically correspond physically to smaller surfactant concentrations, whereas smaller values of * would correspond to larger surfactant concentrations. In principle, the value of * might also be dependent upon air flow rate (dimensionless parameter air ) and upon bubble size (dimensionless parameter 0 ). This could then indicate * being sensitive to the profile of liquid fraction ( vs , which itself depends on air and 0 according to equation (1)). Intuitively a wetter foam (corresponding to higher air and/or smaller 0 ) might be less susceptible to bubble-bubble coalescence than a drier foam would be. However we do not consider such aspects because they were not considered by Keshavarzi et al. (2022b) (only the effect upon * of changing surfactant concentration was explored, as already mentioned).
Clearly * = 1 corresponds to bubble size doubling between = 0 and = base . It is permitted to select foam columns of height very different from base . If we choose a shorter (taller) initial column, there is less (more) variation in bubble size seen than in the base case.
Likewise, over time, the foam column height grows as liquid is removed from beneath the foam, and so we also see more variation in bubble size. However (unlike the work of Keshavarzi et al. (2022b); Grassia (2023)) time variation is not considered here. We look instead just at instantaneous fractionation performance as Hutzler et al. (2013); Tobin et al. (2014) did. A further comment is that the model as written only depends on the ratio * ∕ base and not on the values of * and base separately. Even so, it is convenient to define both parameters, because it means that, when base is set to a specified value (Grassia, 2023) (as will be done here; see Table 1), the value * then gives a more direct measure of relative variation in bubble size for a foam of typical height.

Surfactant flux and effective concentration
In addition to a model for liquid transport we also need a model for transport of surface actives and ultimately enrichment of surface actives. How surface actives behave in a foam, and in particular how they distribute between being dissolved within liquid and being adsorbed on foam film surfaces is a complicated topic (Keshavarzi et al., 2022b), particularly when (Fainerman et al., 2003;Gochev et al., 2013) the surface active material happens to be a bulky molecule like a protein (often targeted in fractionation applications), as opposed to being just a simple surfactant. Here following Grassia (2023) we adopt a highly simplified model: this supposes that if a typical value of surface active concentration within bulk liquid in the foam is given, then a typical value of the amount of adsorbed surface active is also known. Parameter values within the model are informed by Keshavarzi et al. (2022b). Details of this are described in section S2 in supplementary material but we summarise the key points below.
The dimensionless effective concentration ef f is a measure of how enriched the foamate is relative to the feed. It satisfies Here top is dimensionless bubble size at the top and top is liquid fraction at the top. Also Γ 0 * (see section S2) is a dimensionless parameter associated with adsorption of surface actives. Physically Γ 0 * is a ratio between a typical surface concentration and a typical bulk concentration normalised by a length ( ∕( )) 1∕2 , with some geometrical prefactors (see section S2). The value of Γ 0 * is sensitive to foam physical chemistry, but not to parameters like air velocity air , foam column height or bubble size 0 . A typical value is Γ 0 * ≈ 0.025 (see Grassia (2023) and also Table 1) and we use that value throughout.
To give a sense of scale, we mention that one unit of dimensionless concentration in equation (3) might correspond to an actual concentration (Grassia, 2023;Keshavarzi et al., 2022b) of surface active on the order of 0.1 kg m −3 within liquid.
What we anticipate is that increasing foam column height causes ef f to grow, and hence the foam to become enriched in surface actives. This is anticipated to happen for two reasons. Firstly taller columns tend to be drier (Grassia, 2023) ( top is smaller). Secondly taller columns have larger bubbles, and increasing bubble size should cause top to fall even more significantly still. This then is the basis for the expected synergy between foam column height and bubble size acting to enrich the foam. Results are discussed in section 4, but before doing that, in section 3, we consider how the model behaves in a more qualitative fashion.

Analysis of model behaviour
The model to be solved here is defined by equations (1)-(3) along with the boundary conditions already described in section 2.1. The model is easy to solve numerically (see details in section S3 in supplementary material). Before considering numerical model solutions however, one of the ways of analysing the system behaviour (Grassia, 2023) is graphically by sketching, as a function of liquid fraction , various contributions to the liquid flux denoted . It must be remembered however that these contributions necessarily sum to a spatially uniform value thru with liquid fraction at the top then given by top = thru ∕ air . Note that if we set a target top that then determines what the foam column height must be to achieve that particular top value. For the foamate to become enriched, we need top to be significantly smaller than bot , i.e. the foam should become dry at the top. In what follows we analyse a case that is indeed "relatively dry" at the top (in a sense to be made more precise shortly), and then a case which is much drier still.

Case of a relatively dry foam at the top
What we show in the present section is as follows. If there is no spatial variation in bubble size, provided (Hutzler et al., 2013;Tobin et al., 2014;Grassia, 2023) the foam is selected to be relatively dry at the top (in a sense to be specified shortly), then the foam column must become very tall (as will also be explained shortly). On the other hand, if there is spatial variation in bubble size as is the case in the present work, then for the same liquid fraction at the top, the foam column does not need to be nearly as tall.
All this can be explained with reference to Fig. 1. In this figure, the sloping line represents air , the liquid flux convected along with the air flow. Meanwhile the curve (an inverted parabola) represents air − 2 2 , the liquid flux convected by the air and by gravity drainage taken together. The horizontal line thru represents the actual liquid flux carried through the foam, which is spatially uniform. Where the sloping line and the horizontal line intersect defines the liquid fraction at the top of the foam top , where as mentioned top = thru ∕ air . The liquid fraction at the bottom of the foam bot is typically much larger (shown here by a break in the curve).
Meanwhile the difference between the horizontal line and the inverted parabola represents the capillary suction term, which is given by − −1 1∕2 ( 2 )∕ . Suppose now that (as drawn in Fig. 1 for the uppermost inverted parabola) this difference is small, i.e. the peak value peak on the inverted parabola is only slightly less than thru = air top : this then is what we mean by the system being relatively dry at the in a foam column depending on the liquid fraction . The sloping line represents air , the contribution from convection by air. The inverted parabola represents air − 2 2 , the combined contribution from air convection and gravity drainage taken together for a specified value of . This admits a peak value peak at a liquid fraction peak . However the shape of the inverted parabola changes as changes, so strictly speaking the values of peak and peak should change away from their original values (corresponding to an original bubble size 0 ) which are the values marked here. The horizontal line represents the total liquid flux thru delivered through the foam column. The difference between the horizontal line and the inverted parabola represents the effect of capillary suction, which is active when (or more generally 2 ) changes spatially: varies from bot to top , with top corresponding here to the intersection between the sloping line and the horizontal line. As increases, the inverted parabola shifts downwards, so the difference between the horizontal line and the inverted parabola increases.
top. Given ( 2 )∕ is now small (at least for values close to peak , which corresponds to the peak on the inverted parabola), if the spatial variation of is also neglected, then ∕ is necessarily small (again for values close to peak ). We then need a large foam column height for to evolve from bot to top . Indeed much of the height of the foam column is associated with values around peak : as already alluded to, the value of changes only very slowly with there.
Note also that there is a difference in behaviour between a foam which has liquid flux thru greater than the peak peak on the inverted parabola and a foam which has liquid flux less than the peak on the inverted parabola: this is explained in detail in Grassia (2023). In the former case (which is what we consider here) capillary suction is always relevant: as a result the foam has finite height, and also a boundary condition of no slip between liquid and air (as observed by Keshavarzi et al. (2022b)) can be imposed. In the special case when the liquid flux is just slightly above the peak on the inverted parabola, the foam is finite height but tall (Hutzler et al., 2013;Tobin et al., 2014), and we say the foam is relatively dry at the top. On the other hand, when the liquid flux thru is less than the peak on the inverted parabola (not the case we consider here, but a case considered by Grassia et al. (2001); Grassia (2023)), the foam can in principle become arbitrarily tall. High enough up in the foam, capillary suction then eventually becomes unimportant, and the no slip boundary condition cannot then be met.
The situations discussed above, concern spatially uniform values of . The case of spatially varying is different however, because itself grows as evolves. This affects the shape of the inverted parabola. In the figure, it is the uppermost inverted parabola which now corresponds to air − 2 0 2 where (as already mentioned) 0 is the bubble size entering the foam column at the bottom.
There is however now a family of inverted parabolae, one parabola for each . The inverted parabola that is drawn lower down in Fig. 1 corresponds to air − 2 2 with a larger value, remembering that increases here as increases. In effect once has evolved from bot to the original peak say, the liquid flux delivered by air convection and gravity drainage combined has shifted from the upper inverted parabola to an inverted parabola that is lower down. In Fig. 1, this is now quite some distance below the horizontal line thru . As a result, the capillary term need no longer be small, and hence there is no longer a require- Fig. 2. Various contributions to the liquid flux through a foam. In particular the inverted parabola air − 2 2 represents the combined effect of air convection and gravity drainage. The horizontal line represents the total liquid flux through the system thru which is now set at a value below the peak on the inverted parabola. The intersection between the inverted parabola and horizontal line indicates a state with essentially no capillary suction contribution to liquid flux. The value of at this intersection only changes as the value of changes, which in turn changes the shape of the inverted parabola. For large enough , the entire inverted parabola lies below the horizontal line, and capillary terms become relevant again. Liquid fraction varies from bot to top . ment to have slow spatial changes with . In other words there is no longer a requirement for the foam column to become extremely tall for to vary from bot all the way to top .

Case of a very dry foam at the top
In the previous figure (Fig. 1) we considered a thru value slightly greater than peak (the peak on the inverted parabola or, more precisely, slightly greater than the peak on the inverted parabola with = 0 ). Values of thru smaller than this were not considered. Provided varies spatially, there is however now a mechanism to have thru even smaller than the aforementioned peak . This then is what is termed a very dry foam.
First consider (see Fig. 2) the issue with selecting thru below the peak of the inverted parabola but for spatially uniform . At the given thru (horizontal line) the value of evolves from bot to a value at an intersection point between the horizontal line and the inverted parabola (specifically the value at the larger of two intersection points: see section S4 in supplementary material for details).
The value of cannot evolve any further than that, even if the foam column becomes very tall (Grassia, 2023;Grassia et al., 2001). Hence a boundary condition top = thru ∕ air at the intersection between the horizontal line and the sloping line in Fig. 2, cannot be met, at least when is spatially uniform.
In the case with even slowly spatially varying , the situation is different. The value of still becomes held up (i.e. no longer changes much with ) in the neighbourhood of the intersection between the horizontal line and the inverted parabola. However the inverted parabola itself now changes as changes, hence as Fig. 2 shows, the value at the intersection likewise changes (as already mentioned, the resulting formula for is given in section S4 in supplementary material). For each value, we can also identify using equation (2), the value needed to attain that particular . A predicted relationship between and then follows, which is accurate provided is slowly spatially varying, i.e. provided the parameter * is small (again section S4 provides details).
There will be a certain value of (or equivalently a certain value) at which, in Fig. 2, the entirety of the inverted parabola falls below the horizontal line. The value of is no longer held up at the inverted parabola, and so is free now to evolve all the way down to top .
To summarise, provided bubble size varies only slowly with , three regions are expected in the versus profile. In one region (lowest down in ), changes from bot to a value lying on an inverted  We emphasise that the above approximation involves an assumption of small * values (i.e. slowly spatially varying bubble size). Different sets of approximations apply in the limit of larger * values: see sections S5-S6 in supplementary material for details. These latter approximations rely on the observation made earlier that 2 is a dimensionless measure of the Plateau border cross-sectional area. This turns out to be a convenient quantity to analyse, because to an extent, rapid spatial decay in with increasing height is compensated by rapid spatial increase in bubble size . As sections S5-S6 explain, it is then possible to identify which terms in the governing equations are the dominant ones, and approximate accordingly.

Results
This results section presents the predictions of the model used here. It starts off by considering how liquid fraction varies with position within a given foam column (section 4.1). It then moves on to look more globally at how foam column height affects liquid fraction at the top (section 4.2), bubble size at the top (section 4.3), and effective surface active concentration in the foamate (section 4.4).

Profiles of liquid fraction across the foam
Profiles of versus for various top values are shown in Fig. 3. These data are for a fixed bubble size (i.e. * = 0), so correspond to the case already studied in Grassia (2023). However they are included here so as to be able to contrast with variable bubble size cases later on. The main observation is that as top falls, and hence thru = air top likewise falls, the foam column becomes taller. In fact we know (see also Fig. 1) that thru approaches a certain limiting value peak as the foam column becomes very tall indeed. The smallest permitted top (i.e. lower bound for top ) is known (Grassia, 2023) to be air ∕(4 2 ) with = 0 when * = 0. For the particular parameters selected here (see Table 1), this leads to a value just slightly below 0.008.
Note moreover that the domain in Fig. 3 (and some later figures also) is only plotted for ≤ 0.05, although in principle can extend all the way up to bot = 0.36. For greater than about 0.05 though (i.e. near the bottom of the foam column), the various different profiles with different thru all collapse onto the same curve. This is in line with predictions of Grassia (2023), and follows in fact from equation (1). Given that air is typically a small parameter (much smaller than 2 , see e.g. Table 1), then provided is significantly larger than order air ∕ 2 , the dominant balance in that equation is between gravity drainage downwards and capillary suction upwards (Grassia, 2023).  2. An approximate profile is also shown assuming that solutions for lie on an inverted parabola in versus space (see section S4), but the shape of the inverted parabola evolves as bubble size varies with .
A so called equilibrium profile then results (see also equation (S15) in the supplementary material). In effect, for the parameter values selected here, there is no need to examine profiles closely for values greater than about 0.05 as data tend to collapse together as already mentioned. Now we move on to consider cases with nonzero * .
In Fig. 4 some cases with a small but nonzero * are shown, specifically * = 0.2. Again we see that reducing top (hence reducing thru ) causes the foam column height to increase. More importantly though we can now attain top values much smaller than would be permitted with * = 0 and still have a less tall foam column: a profile with * = 0 is also shown in Fig. 4 for contrast.
In section 3.2 we argued that when thru and * are small parameters, the profile of versus can be divided into three regions: a lowermost region, an intermediate region (for which an approximate formula for the profile should be available) and an uppermost region. This scenario is verified in Fig. 5. We present the actual versus profile, and the approximate profile (for details of the approximate profile see section S4). As expected, the actual versus profile approaches the approximate one, stays close to it over just part of the spatial domain (specifically for intermediate values), but then departs again.
Although Fig. 5 considers just one particular nonzero * value, it is possible to argue that the height of the region of intermediate values should scale inversely with * . This follows because the height of this region is governed (as section 3.2 explains) by the entirety of the inverted parabola air − 2 2 becoming smaller than thru , whereas to effect a specified change in the function air − 2 2 , it is necessary to change by a specified target amount. According to equation (2), this then requires to scale inversely with * . Equations (S11)-(S12) Fig. 6. Profiles of liquid fraction versus for * = 2 and various top . Approximations, namely an equilibrium profile and an improvement upon it (discussed in detail in the supplementary material; see section S5), are also shown.
in the supplementary material corroborate this. We have not however presented data with other * values within Fig. 5 because reducing * would take us outside the domain of * contemplated by Keshavarzi et al. (2022b), whereas increasing * would reduce the quality of the fit: we are dealing specifically with a small * approximation here.
The data we have examined to date were obtained for a small * value (bubble size changing only slowly with position). It is also of interest to consider the behaviour with much larger * (bubble size changing more rapidly with position). In Fig. 6 we show profiles of versus for * = 2. Two different values of top are considered. The profiles with different top are actually almost the same over much of the domain (remember the domain extends all the way to bot = 0.36, although we have not plotted all of it here). There is however some sensitivity to top close to the top. It is clear moreover that through having large * , very dry foams (low top ) can now be obtained even over comparatively modest heights (contrast Fig. 6 with Fig. 3 through Fig. 5).
Some approximations (see section S5 for details) are available for the shape of the profile: these are plotted also in Fig. 6. The first approximation follows from the observation, already alluded to earlier, that near the bottom of any foam column (or equivalently over the entirety of the foam column in the case of a column that is not too dry at the top (Grassia, 2023)), the dominant balance is between gravity drainage downwards and capillary suction upwards. As mentioned, a so called equilibrium profile then results. This relies (Grassia, 2023), as alluded to earlier, upon remaining significantly larger than air ∕ 2 , or equivalently upon air ∕ 2 remaining significantly smaller than . Note also that if increases moving up the foam column (as is the case here), then it is more likely that air ∕ 2 will remain small.
An improved approximation can be obtained by retaining the effect of convection by air within the profile: this is also plotted. Details of what this improved approximation involves can be found in section S5. This now includes regions of the profile in which need not necessarily be much larger than air ∕ 2 , although regions right near the top where can become very large indeed (with and top values then being very small) remain excluded. Additional discussion of the data plotted in Fig. 6 can be found in section S7 in supplementary material.

Foam liquid fraction at the top
Thus far we have presented results in the form of versus profiles for various top and * values. Now we turn to consider results more globally by plotting liquid fraction at the top top versus the entire foam column height for various * . Data are presented in Fig. 7.
As can be seen in Fig. 7, making the foam column taller clearly also makes it drier (i.e. increasing , reduces top ). However making bubbles larger is also known in itself to make the foam drier as Grassia  (2023) showed. Hence in a scenario in which foam column height and foam bubble size are coupled (i.e. for nonzero * as we have here), a synergy can emerge: making the foam column taller also makes the bubble size larger, and the combined effect then makes the foam much drier than before. Moreover the stronger the coupling between foam column height and bubble size (i.e. the larger the value of * ), the drier the foam can become (i.e. top reduces further for large * ). A final observation here is the qualitative difference between the * = 0 case and the cases with * ≠ 0. When * = 0 there is a finite lower bound for top even in cases for which becomes arbitrarily large. Keeping * = 0, but increasing the bubble size (by a factor of 1.4, a case also considered by Grassia (2023); in Fig. 7 see the curve labelled "lge ") reduces the value of this lower bound, but a bound remains present. In the case of nonzero * however, there is no such lower bound: making the foam column taller also makes the bubbles larger, which helps with reducing top .
It is clear that increasing * pushes top versus curves to the left in Fig. 7. We can compensate for this by replotting against 1∕2 * . As Fig. 8 shows, this collapses together the data, particularly the data at large * . Data for smaller * however do not collapse quite so well: they tend to be too far to the left for larger top and too far to the right for smaller top . The collapse at large * is however in line with predictions from various approximate formulae in the large * limit: see section S6 in supplementary material.
Recall that (see e.g. Fig. 6) it was possible to generate various approximate profiles for versus . These same approximations can also be used to generate approximate formulae for top versus . The formulae involve either an equilibrium gravity-capillarity balance, or else an improvement upon that taking account of liquid being convected by air (see details in section S5). These approximations (dotted lines) are plotted in Fig. 9 along with the original data (solid lines). Various Clearly at any given * , the approximations perform best when is small (i.e. when top is comparatively large), but less well when increases (i.e. when top is smaller). That said when * is large, the approximate formulae remain valid down to rather smaller top . In fact the approximate analysis in section S5, is more reliable when top remains somewhat larger than air ∕ 2 . This is however relatively easy to achieve when * is large because increases and hence air ∕ 2 decreases.

Bubble size at the top
Thus far we have considered the effect of foam column height on liquid fraction at the top top . However top is not the only physical parameter of interest here. It is also of interest to ask how bubble size at the top top and effective surface active concentration ef f vary: we deal with top variation here and ef f variation in the next section.
In fact the variation of top with is a rather trivial straight line relationship given by equation (2) with = and = top . This is slightly obscured in Fig. 10 by using a logarithmic scale for , but a linear scale for top . That has been done to show that for these data even though varies by up to three orders of magnitude within the figure, the value of top itself only varies by at most an order of magnitude or so. Note also that the curves with higher * within Fig. 10 terminate at lower values than the curves with lower * do. This is because (for nonzero * at least, for which there is no formal lower bound on top ) we have elected to terminate all curves at the same value of top , specifically at top = 0.0005: this value was chosen arbitrarily but is considered to be in the regime of a very dry foam. Curves with higher * reach a given top at a lesser , and that then limits how much top can grow.

Effective concentration
In Fig. 11 we show how effective concentration of surface actives ef f (as given by equation (3)) varies with foam column height . Data for various * values are shown.
The first observation is that the * = 0 case is qualitatively different from the cases with nonzero * . When * = 0, the value of ef f saturates at an upper limit even when foam column height becomes arbitrarily large. This follows from equation (3) remembering that is now fixed at a value 0 . Since (see section 4.1) it is known (Grassia, 2023) that top can fall no lower than air ∕(4 2 ), this limits how much equation (3) allows ef f to grow.
When * is nonzero however, the value of top can fall much lower than this, and so ef f can become very large indeed. This follows because, as already alluded to, there is a synergy between increasing foam   column height and increasing bubble size. This synergy leads to very dry foams (i.e. very low top ) which according to equation (3) are highly enriched.
As already mentioned, we have chosen here to terminate the curves for nonzero * when top = 0.0005 in each case. According to equation (3) this then means we terminate at a ef f that scales inversely with top . Since the cases with larger * terminate at larger top (see Fig. 10), this implies that ef f terminates at a slightly lower value as * increases. Of course this lower terminating value for ef f is associated with a much smaller , i.e. a greatly reduced foam column height. If we compare instead ef f on the various curves in Fig. 11 at a fixed (rather than at a fixed top ), it is then clear that the larger * case leads to higher effective concentration.
Just as we did with Fig. 8 (which collapsed together different top curves) it is possible to collapse different ef f curves together by rescaling the axis. In Fig. 12 we rescale the axis as 1∕3 * , which collapses together data reasonably well, at least for large * . As was the case for Fig. 8 the collapse is less good when * is small. Reasons why collapse is expected here when * is large, and why it involves 1∕3 * (and not 1∕2 * as in Fig. 8) are explained in section S6.

Conclusions
In the context of foam fractionation, we have analysed a simple model for transport of liquid through a foam column. The model predicts the amount of enrichment of a surface active material that thereby results. Within the model, the liquid flux through the foam column is itself spatially uniform, but spatial variations in bubble size are admitted. A parameter * was used to quantify this bubble size variation. In effect * measures the relative change in bubble size from bottom to top in a foam column of typical height.
A qualitative difference was found between systems with no bubble size variation whatsoever ( * identically zero) and systems with slow spatial variation in bubble size ( * small but finite). In the former case there was finite minimum liquid flux that the system could deliver and still meet the boundary condition at the top (no slip between liquid and air). This minimum liquid flux was simultaneously the maximum (the so called peak flux peak ) that could be delivered by air convection and gravity drainage alone without the aid of capillary suction. For small but finite * on the other hand, there was no formal minimum flux that could be delivered, and related with this, no limit upon how dry the foam could become at the top. The aforementioned peak flux peak now in effect decreased moving upwards through the foam column and so, for a tall enough foam column, even a very small imposed liquid flux would eventually exceed peak . A very dry foam leading to a highly enriched foamate could then result although foam column heights needed to scale proportionally to −1 * to achieve that. Turning towards much larger * values, good enrichment of surface actives could be achieved even with comparatively short foam columns. Over much of the foam column, an approximate equilibrium between gravity drainage and capillary suction was found to occur, although the approximation could be improved by incorporating air convection effects. In the limit of large * , the liquid volume fraction at the top top was shown to be a rapidly decreasing function of foam column height . Having small top (i.e. very dry foam) then ensured greatly enriched foamate. Equivalently the foam column height needed to attain a particular top scaled like −1∕2 * in the large * limit. The foam column height needed to attain a given level of enrichment (measured here via an effective concentration ef f of surface actives) scaled meanwhile like −1∕3 * . Overall we identified a synergy between increasing foam height (which even on its own dries out the foam) and increased bubble size at those increased heights (which then leads to much drier foams). Exploiting this synergy is likely to be useful in design and operation of foam fractionation systems, remembering that the drier the foam, the more enriched in surface actives the foamate becomes. It is expected to be particularly relevant in fractionation of comparatively dilute systems for which bubble size might well increase quite significantly moving up through the foam column.
As well as reflecting upon what has been achieved here, it is also worth highlighting what has not been considered. The present work has looked at just instantaneous behaviour of a fractionation column, but previous work demonstrated that fractionation performance evolves over time as the amounts of liquid in a fractionation column falls, and the height of the foam in the column grows to compensate (Grassia, 2023). In view of that, the synergy identified here between increasing foam column height and increasing bubble size is likely to become more significant over time.
P. Grassia and C. Torres-Ulloa Another important point, is that this study has been modelling based: the models still need to be checked against experiment. In particular the models used here for the assumed adsorption behaviour of surface actives have been greatly simplified compared to adsorption behaviour of complex molecules like proteins (Keshavarzi et al., 2022b). In respect of that, it is worth reflecting that the models herein make predictions both of liquid fractions and surface active concentrations/enrichments. Measurement of liquid fractions (typically involving determination of volumes and masses) are arguably easier to achieve than measurement of surface active concentrations (which might involve determination of surface tensions). That said, the models for liquid transport used here are not tied to a particular type of surface active: they could be applied e.g. to a protein (often targeted in flotation applications) but equally they could apply to a simple surfactant. To test the model it may be easier in the first instance to carry out experiments with a simple surfactant, and make measurements of liquid fractions. Once those aspects have been tested experimentally, it should then be possible to progress on to different types of surface actives, including measurements of concentrations/enrichment.

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.

Data availability
Data presented here can be obtained from simple numerical integration of governing equations using parameter values that are provided in the text. Data files are provided with the supplementary material. Data analysed in this work can also be obtained from https:// doi .org /10 .15129 /431c4ad8 -bcf1 -42ec -8fca -e62949b97b76.