Equilibrium state model for surfactants in oils: Colloidal assembly and adsorption

An equilibrium state model addressing the aggregation and adsorption of colloidal assemblies in apolar solvents (oils) via monomer exchange is presented. The model is based on the previously reported step-wise aggregation response of fatty acids and monoglycerides in bio-oils, and captures surface crowding via scaled particle theory (SPT). The sensitivity of key observables - mean aggregation number, adsorbed surfactant amount, and free monomer concentration - to model parameters is demonstrated. Fits to molecular modelling based aggregation and adsorption data of oleic acid and monoolein reveal that the model accurately reproduces chemically speciﬁc aggregate exponential distributions in both bulk and surface phases, even outside of its parameterization conditions. A biased state model, where the initial bulk aggregation step (dimer formation) differs from other steps results in a notable improvement in accuracy. Fits to various phospholipid adsorption isotherms demonstrate the applicability of the model to isotherm type experimental data. The ﬁts reveal either monolayer or aggregate like adsorption structures, depending on surfactant head group charge. The presented model provides an easily accessible, computationally feasible means to estimate colloidal assembly and adsorption in oil environments, and enables assessment of surfactant aggregation propensity and adsorption energetics.

Previously, the kinetics and energetics of adsorption in bio oils on a variety of organo-clay, activated carbon, and mineral surfaces have been investigated using adsorption isotherms measured via quartz crystal micro balance (QCM) [18], interferometric surface force apparatus (SFA) [19,20], rheometry [21], colorimetry [22][23][24][25], and Fourier transform infrared spectrosopy (FT-IR) [26].Furthermore, the structures formed at the oil -solid interface have been characterized using sum frequency spectroscopy (SFS) [27] and X-ray and neutron reflectometry [28][29][30].The findings reveal that adsorption in vegetable oils occurs mainly via the functional groups of chemical species.In addition to adsorbing impurities, the polar ester groups typically found in triglycerides -the major component of most vegetable oils -promote their competitive adsorption on surfaces [18,19,23].Adsorption strength and adsorption capacity is largely determined by hydrogen bonding [18,19,22,31].Adsorption corresponding to monolayer coverage with both perpendicular and parallel molecular orientations depending on molecule chemistry, particularly the sterics of the hydrophobic moieties, have been observed [19,27,30,32].Adsorption is also affected by the effect of moisture and wetting of the surface.Furthermore, the mono-and multilayer adsorption structures observed in dry solvent environments are disturbed by the presence of trace water [33].Additional wetting has also been reported to increase the aggregate structures present at the adsorbent surface [33].
An important contribution to the adsorption of bio-oils arises from intermolecular interactions between adsorbates, particularly in the form of spontaneous self-assembly [18,[45][46][47][48][49].While most bio oils have triglycerides as their majority molecular component, they also contain a diverse range of amphiphilic compounds; the precise composition of which depends on oil feedstock and processing methods [50][51][52][53][54][55].Most naturally present amphiphilic compounds in bio-oils consist of one or more hydrophobic hydrocarbon tails connected by a significantly smaller hydrophilic head group.The head group can be either non-ionic (e.g., monoglycerides) or carry a net charge (e.g., phospholipids).In an apolar solution, these amphiphilic surfactant species self-assemble into reverse micelle-like aggregates, with a hydrophilic core surrounded by a hydrophobic corona [31,[56][57][58][59][60][61][62].The self-assembly of surfactants in apolar solvents differs from that of typical aqueous surfactant solutions: in apolar bio-oil like environments, both monoglycerides and fatty acids exhibit a step-wise aggregation regime which results in an exponential aggregate size distribution [58,31].Additionally, surfactants with significant hydrogen bonding capability, such as monoglycerides or phospholipids, can deviate from an exponential size distribution by showing preference to a given aggregate size regime.Such biased aggregation appears especially at higher surfactant concentrations, where larger aggregates are favoured [58,31].For example, phospholipids -a strongly aggregating surfactant -self-assemble first to cylindrical aggregates and then undergo a further morphological change to spherical reverse micelles with increased phospholipid concentration [63], increased temperature [64,65], or the presence of water or cosurfactant, such as fatty acid [64,45,66].
The aggregation and adsorption response can also be considered by theory and modelling.Ideally, thermodynamic and equilibrium state surfactant aggregation and adsorption models include hydrophobicity driven aggregation, solvent entropy, adsorbentsurfactant chemistry and molecular interactions, as well as the effect of surface heterogeneities (site specific adsorption energies) [67][68][69][70][71].For polar, aqueous surfactant solutions, where physisorption is driven mainly by hydrophobic interaction, the effects of dimerization [72], aggregation [67,[73][74][75][76], differences in surfactant orientation [77], and multi-species adsorption [78,79] on isotherm behaviour have been investigated.To our knowledge, comparable models for apolar solvents do not currently exist.This gap is significant, as readily available prediction and interpretation to adsorption from oils would aid in a number of applications, including oil purification, flotation, and lubrication.
To fill the gap, we propose and implement here an equilibrium state model for describing bulk aggregation and surface adsorption of surfactants in apolar solvent environments.We implement the model, parametrize it against MD-simulated aggregation and adsorption data of mono-and diglycerides from model vegetable oil [31], and assess its robustness and predictive ability.We also test the model performance against experimental data of phospholipid adsorption onto acid activated sepiolite (from rapeseed oil) [14], in particular assessing model sensitivity.For all examined cases, we extract chemically specific model parameters for generalization and comparison.The significance is a wellcharacterized, robust, aggregation and adsorption model for surfactants in apolar solvent environments -to our knowledge, comparable models do not currently exist despite the high technological relevance of both aggregation and adsorption from apolar media.The model is easily generalizable, and provides robust predictions for practical surfactant in oil adsorption systems.

Methods
The thermodynamic equilibrium state model for capturing surfactant aggregation and adsorption in apolar solvent consists of one part capturing bulk solution step-wise growth of aggregates, and another part describing aggregation on a comparatively more crowded adsorbent surface.For simplicity, adsorption is treated as a reversible process and surface heterogeneities including specific adsorption sites are disregarded.The bulk and surface aggregation processes are coupled via monomer exchange between the bulk and surface aggregation models, i.e. all adsorption and desorption takes place as monomers.A conceptual sketch of the reaction equilibria is provided as Fig. 1.The model employs as important parameters the bulk and surface monomer aggregation equilibrium constants K b and K s , and also the free energy of monomer adsorption bDG ads , where b ¼ 1 k b T ; k b is the Boltzmann constant, and T is the absolute temperature.
Following previous findings of stepwise aggregation for surfactants naturally present in plant oils and also cyclohexane type strongly apolar solvents [31,58], the model contains an assumption that the aggregation in bulk solution takes place via reversible stepwise monomer additions.Monomers can also adsorb onto the adsorbent surface, where they are treated as space filling hard discs.On the surface, monomers may further reversibly aggregate to form larger discs.We construct two models: in the simpler model, the same monomer equilibrium bulk aggregation constant K b and surface aggregation constant K s are employed for aggregates of all sizes.However, as the initial aggregation steps are prone to deviate from a stepwise exponential distribution [58], we also consider a model in which the dimerization step in bulk aggregation is governed by its own bulk aggregation equilibrium constant, K b2 , with all subsequent steps governed by K b .This modification accounts for the bias in aggregate size that often arises from increased surfactant hydrogen bonding capability or head group charge [58,80,81], interplay of hydrocarbon tail sterics and head group [82], or the addition of small polar additives, such as water, to the oil [45,57,59,60,64,[82][83][84].For simplicity, surface aggregation is assumed to be less sensitive to the initial aggregation steps, i.e. same surface aggregation equilibrium constant K s is used for all surface aggregation steps.
Aggregation within the bulk solution is described by an analytical solution to the aggregation model proposed by Vierros et al. [58].The equilibrium constant K b for dimer formation 2s 1 s 2 can be expressed as where ½s 1 and ½s 2 are the equilibrium bulk concentrations of monomers and dimers, respectively.Let ½s n then denote the concentration of aggregates with aggregation number n.Hence, ½s 1 ¼ c 1 and ½s n ¼ c n =n, where c n is the concentration of monomers in an aggregate of size n.We generalize the expression of K b for bulk aggregation reaction s 1 þ s nÀ1 s n as Total bulk monomer concentration c bulk is then given by where n max is the maximum aggregate size.The formulation above treats the assembly steps as being thermodynamically equal for all n.Alternatively, as the dimerization equilibrium quite often differs from that of the following monomer addition steps, an explicit equilibrium constant K b2 may be introduced to describe the dimerization reaction 2s 1 s 2 .In this case, K b2 is used merely for the dimerization step and a separate K b is introduced as the equilibrium constant for all consecutive aggregation reaction steps beyond dimer formation, s 1 þ s nÀ1 s n for n > 2. Following the notation of Eq. 2, K b2 and K b can be expressed as This definition of K b2 and K b leads to an overall surfactant concentration Eq. 6 simplifies to Eq. 3 when We assume that the adsorbent surface is more crowded than the bulk solution.Although the concentrations of bulk phase surfactants in apolar solvent environments are generally low, the significantly enhanced concentrations typically occuring at the surface suggests that accurate modelling of the adsorbed phase is important.Here, we assume the adsorbed aggregates to be hard disks, and employ scaled-particle theory to determine thermodynamic properties.Scaled Particle Theory (SPT) [85] has previously been widely used to describe thermodynamic equilibria of adsorption problems, including binary mixture [86], small molecule [87], and protein adsorption [88,89].On the adsorbent surface, aggregates of size n are treated as space occupying hard discs of radius R n ¼ an 1=2 , where a is a fitted constant.In this work, based on fitting to aggregate geometries for palmitic acid and monopalmitin in Ref. [58], we use parameter values a ¼ 0:201 AE 0:008 nm (R 2 ¼ 0:994) for oleic acid and a ¼ 0:353 AE 0:016 nm (R 2 ¼ 0:992) for monoolein.The surface aggregates are treated as a hard disc mixture, with disc-disc interactions defined by a pair potential uðr ij Þ: where r ij is the collision radius of the discs.Here, where R i and R j are the radii of the interacting discs.The discs (aggregates) are not allowed to overlap on the adsorbent surface, which imposes an upper limit on their number density.
SPT is based on the reversible isothermal work WðRÞ needed to expel all discs from a circular cavity of radius R at a given surface occupancy.A Boltzmann factor relates WðRÞ to the probability P 0 ðRÞ of finding a circular cavity of radius R free from any portion of any disc where k B is the Boltzmann constant, and T the absolute temperature.
Our approach is similar to that used by Brusatori and Van Tassel [88] to describe protein adsorption from a binary mixture.However, here we extend the SPT treatment to account for n max differently sized discs and to describing surfactant aggregation and adsorption in oils.To formulate the appropriate SPT approach, we begin with the exact value of WðRÞ for R 6 0: where R n is the radius of the disc corresponding to an aggregate of size n and q n is the surface number density of the discs of that size.
Eq. 9 results from the fact that a cavity with R 6 0 can be blocked by at most one disc.For R > 0, we use a Taylor series expansion for WðRÞ: where W0 is the first order derivative of W with respect to R and p is the hard disk pressure of the 2D disk mixture.To derive p, we follow Brusatori and Van Tassel [88] in noting that the excess chemical Here, where l ig i is the chemical potential of an ideal gas of species i.By employing the Gibbs-Duhem relation to obtain the derivative of p ex and then integrating the resulting equation with respect to q j , the pressure term p is obtained: Now, U n is defined as the probability of finding a cavity of size R n , that is, U n ¼ P 0 ðR n Þ ¼ expðÀbWðR n ÞÞ.Assuming stepwise aggregate growth on the adsorbent surface, in analogy to Eq. 2, the equilibrium constant expression becomes where K s is the equilibrium constant for monomer surface aggregation.
To couple bulk with surface aggregation, a condition on the where bDG ads is the monomer change in free energy of adsorption.Previously, we have calculated comparable free energies of adsorption for fatty acids, monoglycerides, and phospholipids from model vegetable oil solutions using molecular dynamics simulations [31].
Such data can potentially be utilized for model fitting to comparable systems.Combining Eqs. 13 and 14, and limiting the total number of surfactant monomers in the system to N tot to give a closed solution, leads to the following set of non-linear equations: where A is the total adsorbent surface area, and V is the bulk solution volume.
For practical use of the model, we note that the total adsorbent area A and bulk solution volume V, and specificaly the ratio A=V, are system dependent constants.The three main parameters of the model, bDG ads ; K b , and K s , are system independent and can be obtained by, e.g.fitting to existing data points from other, but related systems.Additionally, the aggregate radius R n should be matched with surfactant type, that is, fit to surfactant and aggregate structural data, such as radius of gyration.
The above set of non-linear equations of Eq. 15 is solved for q 1 ; q 2 ; . . .; q n and c 1 using a trust-region-dogleg equation solving algorithm.Penalty functions are applied to imaginary and negative solution roots to steer the solver towards positive and real roots.The concentrations c 2 ; c 3 ; . . .; c n are subsequently calculated based on Eq. 2, or Eqs. 4 and 5, depending on the model.A genetic algorithm was used for parameter fitting and is detailed in Supplementary Information.

Results and discussion
First, we examine the sensitivity of the proposed adsorption model to variation of bDG ads ; K b , or K s .The purpose is to map robustness ranges of the aggregation and adsorption responses for a variety of surfactant species, solvents and adsorbates.Fig. 2a presents the average bulk (solid lines) and surface (dashed lines) aggregation number when each of the three parameters is varied, while keeping the two other parameters constant.Correspondingly, in Fig. 2b, c, and d, the surface concentration of adsorbed surfactant and bulk concentration of total surfactant and free monomer at equilibrium are presented.
Data of Fig. 2a show that the average surface aggregate size responds strongly to K s increasing from zero, i.e. stronger aggregation propensity at the surface.Increasing the adsorption free energy of the monomers bDG ads elicits a similar sigmoidal response in increase surface aggregate size.
The saturation of the surface at high bDG ads or K s is observed in the amount of adsorbed surfactant in Fig. 2b.Surface saturation and extrapolation of the increasing average surface aggregate size curve leads to an adsorption setup corresponding to monolayer coverage on a completely crowded surface.Conversely, an increase in K b results in a depletion of surfactant from the adsorbent surface and a reduction in average surface aggregate size towards only monomers remaining on the adsorbent surface.
For bulk aggregation response, the effect of the three model parameters is, as expected, reversed.In brief, the bulk aggregation response is dominated by K b , with an increase in K b corresponding to the formation of larger bulk aggregates.Now, however, as the bulk solution can be considered to be much more dilute in terms of assemblies than the surface (and indeed no volume exclusion effects are considered for the bulk solution phase here), the bulk solution does not become saturated by the surfactant aggregates.Instead, the assemblies become larger with increasing K b and the bulk monomer concentration diminishes.
The monomer concentration data in Fig. 2d reveals that an increase in bDG ads ; K b , or K s all decrease the free monomer concentration.For K b , this is due to formation of larger aggregates, as we simultaneously observe an increase in overall bulk surfactant concentration (see Fig. 2c).The monomer concentration response to changes in both bDG ads and K s is linked to Eq. 14.Specifically, bDG ads determines the partitioning of monomers between bulk and surface: a more negative bDG ads value results in a smaller monomer bulk concentration as adsorption is more favorable.On the other hand, an increase in K s drives the formation of larger surface aggregates.This results in depletion of monomers at surface.Eq. 14 corrects for this depletion by partitioning more monomer from bulk to surface.Plateaus at large K s or bDG ads values correspond to surface saturation.
Let us next compare the predictions of the thermodynamic model with perhaps the most common theoretical means used to describe monolayer adsorption, the Langmuir isotherm [90,91].The Langmuir model is a useful description for cases where adsorption occurs to well-defined, identical adsorption sites, e.g.chemisorption.The adsorbate is treated as an ideal gas under isothermal conditions, while the adsorbent is assumed to be composed of distinct adsorption sites, each identical [92][93][94].The adsorption phenomenon is a binding equilibrium reaction, where a free adsorbent species, here a surfactant monomer s 1 , binds with a surface adsorption site A to produce a bound adsorbate-absorbent complex A Á s 1 : Let K l be the equilibrium constant for the Langmuir model binding reaction of Eq. 16.Surface coverage h as a function of free surfactant monomer at equilibrium c 1 is then defined as In the model, surface coverage h is one at monolayer saturation.
Here, we define this Langmuir model upper limit based on hexagonal 2D packing of identical monomer sized discs.This leads to a maximum fractional surface coverage of h ¼ p ffiffiffi 3 p =6, independent of monomer disc size.Following the Langmuir model, K l can be calculated directly from bDG ads within the model.The two main assumptions of the Langmuir model are ideal smoothness of the adsorbent surface (all adsorption sites are equal) and the neglect of lateral adsorbate-adsorbate interactions, such as aggregation.An extension of the Langmuir model, the Moreau isotherm [95] includes an adsorbate-adsorbate interaction term I with surface coverage h defined as: where K m be the equilibrium constant for binding at infinite dilution.The adsorbate-adsorbate interaction parameter I describes the attractive (I > 1) or repulsive (0 < I < 1) interaction between adsorbates on the surface.The Moreau isotherm reduces to the Langmuir isotherm when I ¼ 1.
A comparison of surface coverage by the implemented nonbiased aggregation model to the Langmuir and Moreau adsorption isotherms is presented in Fig. 3. To limit bulk aggregation, in the comparison with the Langmuir model, K b ¼ 0:001.The data reveal that an increased surface aggregation propensity, that is, an increase in K s results in a more crowded surface.However, significant surface aggregation is required for the proposed thermodynamic model to predict similar degree of surface coverage as predicted by the Langmuir model.This is expected as Langmuir greatly underestimates the crowding of real surfaces, as it assumes the filling of independent sites.This surface crowding can be accounted to some extent by the Moreau isotherm via introducing a repulsive term between adsorbates, see Eq. 18.Therefore, at K s ¼ 0, where only monomers adsorb, our SPT-based model predicts much lower surface coverage than the Langmuir model, but roughly comparable surface coverage to a Moreau model with a high repulsive parameter.At higher K s , larger aggregates may exist in the SPT-based model which leads to the equilibrium being shifted toward higher coverage.
Next we demonstrate the model parameterization for a chemically specific molecular system, and assess the performance of the model for specific, existing data sets.First, we take MD simulated bulk aggregate and adsorbed aggregate distributions for oleic acid and monoolein in triolein and adsorbed on a quartz or cristobalite surface.The aggregate distributions have been originally published in Ref. [31].The parameters were fitted by the genetic algorithm protocol described in Supplementary Information.Fig. 4 plots the original and fitted bulk and surface aggregate distributions for oleic acid on quartz and cristobalite.The fitted parameters correspond-ing to the graph have been collected to Supplementary Information Table S1.
Altogether, the model predicts accurately the exponential aggregate distributions for both bulk aggregation phase and the aggregation occurring at the adsorbent surface.Additionally, it effectively captures the concentration dependent shift in the aggregate distributions.However, the fitted model appears to underestimate monomer surface concentration.This likely results from underestimation of bDG ads .bDG ads largely determines the partitioning of surfactant between the bulk phase and the surface due to coupling of the bulk and surface aggregation models via monomer chemical potential, see Eq. 14.Additionally, we note that the MD-simulated distributions are subject to simulation time (here 500 ns, i.e. very long) and finite-size effects [31].These likely result in the decrease in fitted bDG ads and K b with concentration for oleic acid on quartz or cristobalite in Table S1 in Supplementary Information.
The aggregates on the adsorbent surface are currently treated as circular discs.In anhydrous conditions, both monoglycerides and fatty acids form filament-like aggregates in oil, that is, one aggregate dimension is significantly elongated in comparison to the other two [31,58,59].Regardless, the assumption of circular aggregates made in this work provides a good fit to fatty acid aggregation data.This observation is unexpected, likely rising from some cancellation of errors, as flexible filaments could be expected to deviate in their adsorption response from circle-like disks significantly.
Fitting of the model to 15 wt-% oleic acid data from the MD simulations enables us to extrapolate the mean aggregate size at also other concentrations, see Fig. 5.The predicted average aggregate size closely matches the data of the MD simulations from Ref. [31].However, the monomer concentrations have mismatch.At sufficiently high concentrations, the mismatch in monomer partitioning leads to over estimation of the mean aggregate size.Altogether, the data of Fig. 5 show that the presented model extrapolates well to concentrations outside its parametrization concentration.
As shown by data in Fig. 4, the model that has a single bulk aggregation equilibrium constant (see Eqs. 2 and 3) produces a good fit to fatty acid data.Additionally, the assumptions of i) adsorption occurring as monomers, and ii) identical free energies of adsorption for free and aggregate monomers, lead to good agreement with data from apolar solvent solutions of surfactants with Fig. 3. Adsorption isotherms based on the non-biased adsorption model (solid lines) compared to prediction according to Langmuir and Moreau adsorption theory (dashed lines).The data is for K b ¼ 0:001 and bDG ads ¼ À5:0.Surface area is 500 nm 2 and Rn ¼ 0:5n 1=2 nm.MD simulations data originates from Ref. [31].Fig. 5.The average bulk aggregate size predicted by the fitted thermodynamic model with a comparison to average aggregate sizes from MD simulations data originating from Ref. [31].The MD simulations data rise from a variety of concentrations ranging between 5 wt-% to 25 wt-% initial bulk concentration of oleic acid while the fitted thermodynamic model prediction is based on the model fit to the 15 wt-% oleic acid MD data.
relatively small polar head group and limited hydrogen bonding capability, i.e. fatty acids.However, for surfactants with larger hydrophilic headgroups (and increased hydrogen bonding capability) or surfactants with charged groups (e.g.fatty acid soaps or phospholipids), an aggregation bias towards specific aggregate size exists [31,58,81,96].Similar effect arises from water dispersed in bio oils as surfactant headgroup hydration leads to formation of water bridges and aqueous cores in surfactant reverse micelles and other aggregates [45,57,59,61,[97][98][99][100].Although a simplification, the biased aggregation model where dimer formation is differentiated from other bulk aggregation steps, as defined by Eqs. 4 and 5, captures the essence of this bias.
The monoolein data set of Ref. [31] is an example of such enhanced aggregation propensity due to increased surfactant head group hydrogen bonding capability.This also leads to a stronger binding to the adsorbent surface in comparison to oleic acid.Indeed, this is reflected in larger K b and bDG ads values, see Table S1 in Supplementary Information, in the fit to MD data.Additionally, at higher monoglyceride concentrations, larger surfactant aggregates become thermodynamically more favoured compared to initial aggregation steps [58,31,59].In the thermodynamic model parameter fit (Table S1), this shows as a strong concentration dependency of K b for monoolein.
To account for aggregation energy differences between initial aggregation steps (dimerization) and later aggregation steps, we include also an explicit dimerization constant K b2 that may differ from the bulk aggregation equilibrium constant K b for later aggregation steps (n > 2).The equilibrium constants K b2 and K b are defined in Eqs. 4 and 5. Let us next fit this biased aggregation model to the monoolein aggregate distributions.The fitted K b2 ; K b ; K s , and bDG ads values are presented in Supplementary Information in Table S2.Fig. 6 plots the fit of the biased model to the MD simulations data for 10 wt.% concentration of monoglyceride in triolein solvent and adsorbing on quartz.Notably, the fitting results in a lower K b2 compared to K b .This means that the initial dimerization is significantly less thermodynamically favoured than subsequent growth steps.A similar bias has previously been noted by Vierros et al. [58].In the current fit, each point in the MD-simulated aggregate distribution is weighted equally.For a better fit, uneven weights could be applied to better account for the scatter at the large aggregate end of the distribution, see especially Fig. 6a); in the current fit, this scatter in original data set leads to a significant underestimation of both bulk monomer concentration and small surface aggregate concentrations.
Indeed, individually tuned weights for different aggregation steps would bring access to additional resolution in terms of chemical specificity of the modelled system.However, this route is at the cost of simplicity and leads to more complex, parameterization data heavy solutions.Additionally, creating such complex models does not guarantee improved interpretability.
To further generalize the model for surfactants with larger polar headgroups, we next demonstrate the performance and fitting of the model to phospholipid adsorption isotherm data for DOPC, DOPA-HNa, and DOPA-H 2 in rapeseed oil on acid activated sepiolite.The isotherm data measured by inductively coupled plasma spectrometer were originally published by Laatikainen et al. in Ref. [14].The residual water content of 0:2 wt-% in the measured system is expected to promote aggregation of the phospholipids when bound to the polar headgroups although the level of hydration depends on phospholipid type.At low phospholipid concentrations ($ 0:5 mmol/kg), that are well below the experimentally observed CMCs [101], water content between 0:03 and 1:00 wt-% was reported to have minimal effect on overall adsorption [14].While the existence of CMC for surfactants in apolar solvents is debatable [31,58,102], it is well known that water, as well as phospholipid composition and solvent medium influence the aggregation propensity and reverse micelle morphology of phospholipids [45,64,66].Here, step-wise growth of bulk phospholipid aggregates is assumed with no present CMC.We note that while evidence of step-wise aggregate growth exists for weakly polar surfactants such as fatty acids and monoglycerides [31,58,102], evidence of similar aggregate growth for charged surfactants, such as phospholipids, remains sparse in current literature [103].Due to both low phospholipid concentration and low hydration of the phospholipid heads, we choose the simplified model with one bulk aggregation equilibrium constant as the comparison model.For the fit, a surface area of 176 m 2 /g for the acid activated sepiolite is assumed [14].In addition to DG ads and K b , a radius r to describe average molecular area is also fitted.As phospholipids are expected to adsorb as monolayers at hydrophilic solid -liquid interfaces in apolar solvents [14,104], no aggregation at adsorbent surface is included in the fitted model, i.e.K s is set to zero.
The original data by Laatikainen et al. and our model fits are presented in Fig. 7a.The fitted values for DG ads ; K b , and r are as follows: bDG ads ¼ À13:31 (' À37:43 kJ/mol), K b ¼ 9:68, and r ¼ 0:47 nm for DOPC, bDG ads ¼ À11:13 (' À31:28 kJ/mol), K b ¼ 6:62, and r ¼ 0:60 nm for DOPA-HNa, and bDG ads ¼ À9:97 (' À28:04 kJ/mol), K b ¼ 13:42, and r ¼ 0:39 nm for DOPA-H 2 .Notably, the bDG ads values resulting from our model fit are $ 4 -11 kJ/mol higher than those reported by Laatikainen et al. based on a Moreau isotherm fit to the adsorption isotherm data [14].On the other hand, for DOPC and DOPA-HNa, the fitted r values fall between those reported by Laatikainen et al. [14] based on observed surface saturation and measured sepiolite surface area, and molecular surface areas based on packing of the respective phospholipids in a bilayer (0:46 nm for DOPC and 0:50 nm for DOPA-HNa) [105,106].For DOPA-H 2 , the fitted r value falls below that expected for bilayerlike packing (0:49 nm) [106].We conclude from this that the packing of DOPC at saturation is close to that of a bilayer.On the other hand, the sparser packing of DOPA-HNa on the surface corresponds with the expected electrostatic repulsion of the charged head groups.Finally, for DOPA-H 2 , the tighter packing suggest the presence of surface aggregates, even though surface aggregation was ignored in the current fit.
When no surface aggregation is included in the model, the choice of K b has minimal effect on the amount of adsorbed surfactant at equilibrium.Instead, the partitioning of surfactant between bulk and surface is determined by bDG ads while r determines the saturation concentration of the isotherm.The insensitivity of the isotherm is demonstrated in Fig. 7b, where the original fitted K b value for DOPC is scaled by a factor of 0:1; 2, and 10 without significant changes to the predicted isotherm.The significance of this is that if the proposed model is fit to isotherm data, K b can be considered as a free parameter.
The three main model paramters -bDG ads ; K b , and K s -can be directy mapped to experimentally observable properties, specifically the free energy of adsorption and the free energy of micelization (or dimerization).The mismatch between experimentally determined free energies of adsorption -particularly when derived via isotherm fitting -and the model prediction can be significant as demonstrated here.Alternatively, molecular modellling techniques can be used to determine the adsorption energies of single surfactants [31,34,58,107,108] and free energies of aggregation [58,[109][110][111][112].Experimentally determined free energies of micellization [113][114][115] could also be mapped to K b equilibrium constant.The value of K s , that is, the aggregation propensity at the surface, is related to K b and surfactant -surfactant interaction strength.Additionally, it is modified by the effects of surface diffusion and adsorption geometry of single surfactants and preaggregates.
While the model accurately replicates aggregation and adsorption regimes in the examined simple model systems, improvement prospects remain.First, while native concentrations of amphiphilic impurities in typical bio-oils are low, deviation from ideal mixture behaviour should be expected, especially for strongly hydrogen bonding or charged species, or in the presence of water.Such non-ideality could be accounted for by introducing activity-based correction factors to the model.Additionally, the SPT potential could be augmented by an additional Yukawa potential based term to account for electrostatic repulsion, which is relevant for adsorption of charged species, such as phospholipids or soaps.The adsorption of charged surfactant species is also sensitive to the presence of water [14,116], ions [14,117], and charged adsorbent surface moieties [34,101], which all contribute to the screening of electrostatic interactions and therefore packing of the surfactant at the interface.Additionally, surface morphology guides the packing [118,119].Also the treatment of surface aggregates could be improved to better describe e.g., the rod-or filament-like aggregates [31,[58][59][60] formed by species such as monoglycerides and fatty acids in anhydrous apolar solvents.We note, however, that despite the current model assuming circular disks, a relatively good fit is obtained for also such filament-like systems.Indeed, the SPT-based adsorption has been priorly modified for adsorption of particles of arbitrary shape [120].However, the treatment of surface aggregates as two-dimensional space filling discs abstracts the three-dimensional packing of the surfactants.Particularly, the surfactant tail conformation changes upon adsorption can be significant, leading to an entropic contribution [119,[121][122][123].This is now omitted in the model.Finally, the adsorbent surface is assumed perfectly planar here and no information on surface roughness or preferential adsorption sites or conditions is included in the model.Advancements of the model could include considering the roughness and its effect on the adsorption for better match with real, non-ideal adsorbents.Another potential future extension of the model includes implementing a competing adsorbate species into the model, which would most likely capture better the behaviour of real oil systems, as well as provide insight into the diverse equilibrium phases based on e.g., the relative adsorption energies of the two competing adsorbates.

Conclusions
We presented and implemented an equilibrium state model for describing step-wise aggregation and adsorption of surfactants in apolar solvent environments.The model encompasses step-wise bulk aggregation coupled with step-wise surface aggregation via monomer exchange between the bulk and the surface populations.Crowding of the surface is accounted for by an SPT approach.To account for a bias toward specific aggregate sizes in bulk or at the surface due to e.g., ionic groups, significant hydrogen bonding capability, head group hydration, or small polar additives in the apolar solvent environment, we also introduced as the simplest solution an explicit dimerization constant that may differ from the bulk aggregation constant used for the subsequent aggregation steps.We demonstrated the sensitivity and robustness of the models against variation of their parameters, and mapped the performance against both microscopically detailed simulations and experimental adsorption isotherm data of typical surfactants in bio-oils.The presented particle-level thermodynamic approach to describing colloidal assemblies and their adsorption response in low dielectric environments is, to our knowledge, a currently non-existing approach to predict equilibrium response of apolar solvent surfactant solutions.Most importantly, the approach overcomes the challenges associated with molecular level descriptions of such systems [124][125][126].
In summary, we presented here a model capable of predicting practical aggregation and adsorption behaviour in a set of technologically relevant bio-oil surfactant adsorption systems.The model can also be used to assess the likely adsorption structures based on the parameter values relating to surfactant chemistry, aggregation propensity, and adsorption strength.These parameters can be tuned to match specific chemical systems and experimentally accessible observables.For example, packing density on adsorbent surface, mean aggregate size, and variance can be extracted, as demonstrated by the fits to phospholipid adsorption data.Understanding of such surfactant equilibrium assemblies plays an important role in e.g., optimising sorption-based extraction and purification processes.Furthermore, bulk phase equilibrium aggregates may be linked to e.g., rheological properties of microemulsion systems.In total, the presented model brings access to equilibrium predictions and interpretation of the response of a wide variety of surfactant -apolar solvent -adsorbent systems.

Fig. 1 .
Fig. 1.Conceptual scheme of the aggregation and adsorption reaction equilibria considered by the state model.Parameters include bulk and surface monomer aggregation equilibrium coefficients K b and Ks, and monomer free energy of adsorption DG ads .Aggregation may be biased by separating K b into a dimerization equilibrium constant K b2 differing in value from the equilibrium constant K b corresponding to all assembly steps beyond dimerization, i.e. for aggregation number n > 2.

Fig. 2 .
Fig. 2. Effect of bulk and surface aggregation equilibrium constants K b and Ks, and monomer adsorption energy bDG ads on modelled a) equilibrium aggregate size in bulk solution (solid lines) and on adsorbent surface (dashed lines), b) concentration of adsorbed surfactant, c) concentration of surfactant and d) free monomer in bulk phase at equilibrium.Dependency of K b is mapped with keeping bDG ads ¼ À2:0 and Ks ¼ 1:0 constant, Ks with keeping bDG ads ¼ À2:0 and K b ¼ 1:0 constant, and bDG ads with keeping K b ¼ 1:0 and Ks ¼ 1:0 constant.Initial bulk surfactant concentration is 0:5 M, surface area 100 nm 2 , and Rn ¼ 0:5ðnÞ 0:5 .

Fig. 4 .
Fig. 4. The MD simulated (filled symbols) and fitted model (open symbols) aggregate distributions for oleic acid in triolein.Panel (a) presents the bulk and panel (b) the surface aggregate distributions on quartz.Corresponding distributions on cristobalite are presented in panel (c) for bulk and in panel (d) for surface aggregate distributions.MD simulations data originates from Ref. [31].

Fig. 6 .
Fig. 6.Aggregate size distribution (a) in bulk solution and (b) on surface for 10 wt-% monoolein in triolein on quartz as predicted by MD-simulations and fitted model with explicit K b2 and K b .

Fig. 7 .
Fig. 7. (a) Adsorption isotherms (bulk equilibrium concentration vs. adsorbed surfactant per kg of clay) for DOPC, DOPA-HNa and DOPA-H 2 from rapeseed oil onto acidactivated sepiolite at T ¼ 65 C, and water content of 0:2 wt-%.Filled symbols denote the experimental isotherm data originally published by Laatikainen et al. [14], solid lines show the fitted model isotherm prediction.(b) Fitted adsorption isotherm for DOPC plotted using different values of K b .