Modeling surfactant and drop size dynamics in polydisperse liquid-liquid systems with population balances

(cid:1) The high order moment conserving method of classes for solution of population balances is extended to include surfactant dynamics. (cid:1) The proposed method is fast, robust, and very accurate even with a limited number of variables. (cid:1) The novel method is tested with realistic liquid–liquid dispersions with relevant physical closure models for surfactant adsorption, mass transfer and droplet breakage and coalescence rates.


Introduction
Liquid-liquid dispersions are encountered in many Chemical Engineering applications, such as in liquid-liquid extraction, in some polymerizations and other chemical reactions, pharmaceuticals, food industry, oil refining, and in general often when organic and aqueous phases are present simultaneously.Often drop size distributions are of high importance in these systems, either providing surface area for mass transfer or suitable structure in emulsions.Drop sizes and drop dynamics also determine operation of liquid-liquid coalescers (Seader et al., 2011;Paul et al., 2003).Surfactant adsorption to liquid-liquid interfaces is an important phenomenon in many of these applications.Adsorption can be either intentional, when surfactants are added to the system in order to reduce drop sizes and/or stabilize emulsions, or unintentional when surfactants are merely unavoidable components in the processed streams (Tardoz, 2015;Abbott, 2017).
The surfactant adsorption is a dynamic phenomenon covered by mass transfer to and from the interface, equilibrium composition at the surface (typically described by an adsorption isotherm) and dynamics related to the formation of liquid-liquid interface by breakage and coalescence of droplets.Proper models are needed to account for all these phenomena (Chang and Franses, 1995;Maindarkar et al., 2015).For a detailed discussion about surfactant adsorption dynamics to interfaces, see e.g.Ravera et al. (2000) and Miller et al. (2017).
There is rich literature available for detailed simulations of liquid-liquid interfacial behavior at small scales, even in the presence of surfactants (James and Lowengrub, 2004;Muradoglu and Tryggvason, 2008;Xu et al., 2006;Dziuk and Elliott, 2013).Although these surface capturing fluid dynamics approaches can give valuable insight to drop scale dynamics, such simulations cannot in practice be connected to modeling on equipment scale, where numerically resolving each drop surface is simply impossible.However, well-established population balance framework is amenable also to process-level approach, and has been applied to

Contents lists available at ScienceDirect
Chemical Engineering Science j o u r n a l h o m e p a g e : w w w .e l s e v i e r .c o m / l o c a t e / c e s liquid-liquid systems for quite some time.There has been extensive development of drop breakage and coalescence rate models suitable for population balances (Ramkrishna, 2000;Coulaloglou and Tavlarides, 1977;Alopaeus et al., 1999b).Proper combination of population balance methods with detailed numerical and experimental studies on drop size level is probably the most powerful tool for rigorous modeling of process equipment containing two immiscible phases.This can be done with a simplified flow models for the equipment together with material and energy balances, or by applying rigorous CFD tools to model the flow field and its details required for the drop breakage and coalescence rate models.There are also successful efforts to implement surfactant adsorption in these models (Chen et al., 2019;Håkansson et al., 2009Håkansson et al., , 2013)).
In this contribution, a high order moment conserving method of classes (HMMC) developed earlier for population balances and polydisperse drop phase concentrations (Alopaeus et al., 2006b(Alopaeus et al., , 2007(Alopaeus et al., , 2008;;Buffo and Alopaeus, 2017) is extended to capture drop size dependent surface phase material balance for the surfactants.This allows fast and accurate solution method for drop size, bulk material, and surfactant balances, paving a way further to truly designing practical equipment using population balances.

Population balance model for droplet sizes, concentrations, and surfactants
The proposed method is based on the high order moment conserving method of classes (HMMC).The original method was first developed for bubble and droplet population balances with breakage and coalescence, then extended to growth and nucleation, and finally a consistent distribution reconstruction method was proposed (Alopaeus et al., 2006b(Alopaeus et al., , 2007(Alopaeus et al., , 2008)).Later the method was extended to multi-dimensional population balance modeling (Buffo and Alopaeus, 2016), applied to liquid-liquid extraction columns (Buffo et al., 2016), cellulose and hemicellulose chain length distributions (Visuri et al., 2012;Ahmad et al., 2015) and flotation (Basavarajappa et al., 2017).
A fundamental advantage in the high order moment conserving method of classes is that it predicts distribution integral properties (distribution moments) with considerably smaller number of pop-ulation balance related additional variables compared to traditional method of classes, but is able to predict distribution shapes in more details than traditional (and quadrature) moment methods.The method does not require any iterative solutions for solving population balances, enables faster model specific time integration in some relevant cases than with standard solvers (Jama et al., 2020), and is suitable to be implemented in computational fluid dynamic solvers (Jama et al., 2021).Due to this, HMMC can be implemented according to the requirements of each modeling case: with a very small number of categories the distribution moments are already predicted with reasonable accuracy, but if the precise distribution shape is important, the same method can predict it practically without any numerical error in the moments by adding a few more size categories.
The approach presented in this contribution extends simplified pseudo-2D population balance (Buffo and Alopaeus, 2017) for drop size and concentration distribution to include drop size dependent material balances also for separate surface phase.The simplified pseudo-2D method was shown to perform almost indistinguishably as compared to a more rigorous bivariate method for practical cases with mass transfer and reaction, including feed and outlet source terms mimicking an extraction stage of an extraction column.The present further development allows tracking also surfactant dynamics due to mass transfer, drop breakage, and coalescence.It opens new rich research avenues related to multiphase system analysis containing surfactants.This is possible due to fast and simple numerical solution tool.Although the method is presented here primarily for liquid-liquid systems (to avoid cumbersome referencing to all different phase combinations), it can readily be extended to gas-liquid systems and systems with multiple dispersed phases.
In this chapter, the proposed method is explained and described first in its generic form, where closure models related to mass transfer, adsorption and phase equilibrium as well as drop breakage and coalescence models are left undefined, and later in the Example -chapter some of the most typical closures are presented in the context of relevant test cases.Feed and discharge of droplets and corresponding material is not included in these model equations to avoid unnecessary complexity.However, source terms crossing the control volume boundaries (whether differential or The underlying principle of the present distributed material balance model is described in Fig. 1.For the continuous phase, material balances for each chemical component is formulated normally.These material balances can be localized, i.e. written for a computational cell in CFD, or length coordinate in tubular reactor model, or they can be global as in a CSTR model.In each case the system can be assumed to be in steady state or time dependent.Drop sizes are modeled by discretizing the drop size distribution into predetermined size categories, where (local) number concentration of each droplet size is followed.For the dispersed phase (inside droplets), material balances are formulated separately for each chemical component and drop size.Different droplet sizes can thus have different compositions resulting from differences in mass transfer rates and droplet phase mixing processes related to breakage and coalescence.
In this work, similar material balance distribution is here applied for the interface as for the drop phase, and relevant equations presented.It is expected that those components exhibiting surface activity are concentrated at this thin interface.The present model thus follows amount of moles at the interface, calculates surface concentrations based on amount of moles and available surface areas, and assumes that there are two mass transfer films on both sides of the interface.Interface adsorption equilibrium sets concentrations in the immediate vicinity of the interface in both liquid phases (continuous and droplet).Driving force for mass transfer is the difference between bulk phase concentrations and concentration in the immediate vicinity of the interface.Interface itself is here considered to be very thin, so that there is no additional mass transfer resistance in it, but the liquid phases in the immediate vicinity of it on both sides of it are also in phase equilibrium.Schematic illustration describing mass transfer and relevant phase equilibria is shown in Fig. 2, and variable treatment in the discretized balances are illustrated in Fig. 3.
The model equations for population balance, droplet phase material balance, and interface material balance are shown below.
Population balances: Droplet phase material balances: Interface material balances: Continuous phase material balances: There are different ways to model mass transfer fluxes, starting from a simple approach where individual component mass transfer coefficients are multiplied with corresponding driving forces, up to rigorous multicomponent approaches where diffusional interactions are taken into account in the mass transfer coefficients, thermodynamic corrections in the driving forces, electroneutrality requirement in the fluxes as well as convective flux (Taylor and Krishna, 1993).The most important thing in the spirit of the present interface model method is that phase and surfactant adsorption equilibria should be consistently defined for the two bulk phases and the interface.Note analogy in the treatment of population and the two material balances: the discretization matrices A, B and G are the same for all balances, and can be calculated prior to the simulation.These depend on the selected discretization and method order (which moments are selected to be conserved).For details related to construction of the discretization matrices, see the original papers describing HMMC method (cited above).Treatment of droplet and interface material balances are completely analogous for population balances, while mass transfer needs to be defined consistently according to the defined mass transfer directions.All the calculations are fully explicit from the population balance viewpoint, although mass transfer fluxes may require iterative solution  if formulated rigorously (Taylor and Krishna, 1993;Alopaeus and Aittamaa, 2000;Alopaeus et al., 1999a;Alopaeus, 2002).In practical implementation, all the balances can be solved in parallel (within the same for-loops).The material balances are formulated here based on extensive variables (moles), although they could be formulated for concentrations in an analogical manner.In the present formulation, volume changes due to excess volume in mixing, reactions, or other time dependent reasons except mass transfer are neglected.If these turn out to be important, they can nevertheless be added to the present model.The method itself is conserving droplet volume, droplet phase, and surface phase mass exactly: a mass source term in one category or phase is always precisely the same mass source for the others when droplets are breaking or coalescing.This leads to overall mass conservation in the numerical method to be limited only by the time integration tolerances, or finally by floating point calculation errors (''machine epsilon").When the model is implemented in any simulation environment, it is important to verify the implementation with extensive testing to guarantee this.A standard method for this is to compare initial dispersed phase volumes and amounts of material fed to the system to those after the simulations and confirm that no changes are occurring in a steady state with simulations significantly longer than any initial transients.

Closure models
Although the closure models for mass transfer, breakage and coalescence, and other relevant rates are always necessary for the proposed population balance model to work, selection of proper closures is always case-dependent.Some of them, such as thermodynamic basis for phase equilibria and surface adsorption, are very generic, and some of them, such as the applied mass transfer coefficient correlation, may be adopted more freely based on suitable empirical correlations for each system.The proposed modeling framework for population and drop size dependent materialand surfactant balances is very flexible in terms of selected closure models.The model does not require case-dependent derivation of implemented equations for each closure, but instead normal rate and equilibrium models can be used as such.The selected closure models for further numerical illustrations are described next.

Phase and surface equilibria
Mass transfer between surface and continuous phase as well as surface and dispersed liquid phase are essential parts of the model.For surfactants, mass transfer is always calculated with interface composition, although for transfer of components which are not surface active, this is not necessary and normal two-film theory can be used.Surface equilibrium compositions are calculated using Langmuir adsorption isotherm, which is quite often used equilibrium model especially for non-ionic surfactants (Chang and Franses, 1995;Abbott, 2017).In the Langmuir isotherm, there are two adjustable parameters, K L describing surface activity and C max describing maximum surface coverage: Typically, maximum surface coverage does not vary much between surfactants, but surface activity may vary by several orders of magnitude.Langmuir isotherm can be extended to multiple surfactants relatively easily if surface coverages for different components are the same, but in the present example we focus on single surfactant systems.Other adsorption isotherms can be used as well in the proposed method, but as the present work focuses on method development rather than detailed analysis of specific surfactant containing dispersions, this is not exemplified in the further illustrations.
Significant fraction of literature concerning surfactants is focusing on gas-liquid interfaces.For liquid-liquid interfaces, the surfactant adsorption from both sides of the interface follow similar physical steps than in a single liquid case, although mass transfer rates naturally depend on whether the bulk phase is continuous or dispersed.One important aspect constraint is that the adsorption and phase equilibria need to be consistent.Phase equilibria between bulk phases is typically expressed with distribution coefficients (Seader et al., 2011): where the latter equality describes liquid-liquid equilibrium from thermodynamic point of view, originating from equal fugacities of each chemical component in both liquid phases.Here superscripts I and II refer to the two liquid phases.It is customary to set lighter (smaller density) as I and heavier as II, although this is not always a straightforward definition as one phase may be lighter at certain composition range and heavier at other.In all cases, it is important to specify the phases consistently.Quite often one phase is aqueous and another organic.Whichever of these is lighter or heavier on one hand, or dispersed or continuous on the other, needs just to be treated systematically.
In terms of concentrations (as concentrations are typical variables in adsorption isotherms), the bulk equilibrium equation can be written as where c tot refers to the total concentration (molar density) of a phase.If molar densities are included in the distribution coefficient, the equation can be written as There are thus three equilibrium conditions: equilibrium between dispersed phase and interface, equilibrium between continuous phase and interface, and equilibrium between bulk phases.For each component, two of these can be considered independent, and the third one follows from those two.
In many cases surfactants have notable solubility only in one of the bulk phases due to a very high activity coefficient in the other.In such cases bulk equilibria is not relevant from the modeling point of view, and adsorption equilibrium from one phase to the interface is sufficient to model interfacial chemistry.The present model is nevertheless formulated in a generic form without a need for further assumptions regarding solubilities.

Mass transfer rates
Mass transfer rate consists of diffusive and convective contributions.For inter-phase mass transfer, convective contribution and diffusional interactions are often neglected, although for cases where significant amounts of several components transfer between the phases, those should be included (Taylor and Krishna, 1993).When surfactant is the only transferring component, the system can be typically considered dilute (with respect to the transferring components), and a simpler mass transfer model suffices.In the present examples, mass transfer is modeled with a simple mass transfer coefficient -based model as where mass transfer is defined positive for fluxes from continuous to dispersed phase.Mass transfer coefficients are defined separately for each transferring component and drop size, although in the following examples only one component is transferring.Concentrations for each drop size are defined as discussed earlier.Interface concentrations in this formulation refers to such concentration which is in equilibrium with the interface in the immediate vicinity of it.Final adsorption dynamics, i.e.rate of adsorbing surfactant molecules to attach and orient at the interface, is assumed to be instantaneous.
Mass transfer coefficient for the dispersed phase is calculated from a rational approximation of diffusion within a sphere (Alopaeus, 2000), where transient time for diffusion is calculated from the drop breakage frequency as in Laakkonen et al., 2006a, 2006b, 2007and Alopaeus, 2014.For the continuous phase, penetration theory is applied: where terminal velocity of a droplet is calculated with the following correlation (Vignes, 1965;Buffo and Alopaeus, 2017):

Breakage and coalescence rates
Breakage and coalescence rate models for the liquid-liquid dispersion are taken from Alopaeus et al. (2002) with slight modifications explained below.The breakage rate is given by g Note that in the present formulation, interfacial tension is drop size dependent.
For the daughter size distribution, the beta distribution is chosen due to its consistency with the breakage parameters optimized in Alopaeus et al. (2002) Coalescence rate consists of two parts, coalescence frequency and coalescence efficiency.The coalescence frequency originates from Coulaloglou and Tavlarides (1977): Collision efficiency was developed earlier (Alopaeus et al., 1999b) as As can be seen, in the original formulation, impeller power number and diameters are included in the collision efficiency function to better represent the fluctuations in impeller driven flow.In order to generalize the collision efficiency function further, typical values for impeller power number and diameter in laboratory scale stirred tanks (where the drop size distributions were originally measured) are inserted.If values Np = 5 and D i = 0.1 m are used, the following equation results in: where C 0 11 = 2.71446[m À2/3 ].As the original parameter C 11 was not fitted against experimental data, but an order of magnitude was estimated instead, the present generalization is probably not a major source of discrepancy when predicting drop size distributions in various turbulent systems, especially if impeller type and size are relatively close to the original setting used in the parameter optimization.A drawback in this approach is, that the new parameter C 0 11 is dimensional, implying that upon stirred tank scale-up, droplet coalescence efficiency may not depend only on turbulent dissipation rate, physical properties and drop sizes.
Coalescence efficiency seems to be the most controversial part in the drop breakage and coalescence rate models.The efficiency adopted above does not depend on interfacial tension, while the coalescence efficiency adopted by Coulaloglou and Tavlarides (1977) does not depend on dispersed phase properties.Further, the model by Coulaloglou and Tavlarides with their original parameters predicts almost 100% efficiencies for drop coalescences in many cases where droplets are small, characteristic to stirred tank operations.This implies that the optimized coalescence efficiency parameter may not be truly identified during optimization, and the whole efficiency contribution may remain redundant.Recently there have been detailed simulations predicting droplet coalescence phenomena on drop surface deformation level, hopefully leading to practical efficiency correlations which can be used in the population balance models without additional computational burden.However, the physics is very complex consisting of different interactions between the drops (Ozan, 2021;Liao and Lucas, 2010).The proposed population balance solution framework could help to improve these models, as drop size dependent surfactant properties are available during the model solution for further closure model development and parameter optimization against available experimental data.
As the adopted drop coalescence efficiency does not account for effect of surfactants, the correlation proposed by Håkansson et al. (2009Håkansson et al. ( , 2013) ) is used to consider effect of coalescence inhibition by the presence of surfactants: In the present example, coalescence efficiency is multiplied with this inhibitory factor.
Prescence of surfactants also affects interfacial tension, which is an important parameter especially in the drop breakage functions, where interfacial tension is typically the most important factor stabilizing the droplets against breakage caused by turbulence.Interfacial tension can be calculated from Gibbs' adsorption equation, which reads for Langmuir isotherm and non-ionic surface concentration as This is one form of the Szysztowski equation.Additional constraint for Langmuir isotherm is taken here as the critical micelle concentration (CMC).It is assumed here that all surfactant molecules in the bulk phase above this concentration form micelles, and consequently, surface concentration cannot exceed a threshold value which is determined by CMC and corresponding interfacial tension.This assumption is used also for the mass transfer modeling, where driving force for mass transfer is limited with bulk CMC.Interfacial tension -actual concentration dependency is illustrated qualitatively in Fig. 4.
It should be noted here that the present assumption where CMC is limiting surface coverage, leads together with the assumed coalescence rate suppression model, that even with surfactant bulk concentrations well above CMC, some coalescence will take place although coalescence rate is greatly suppressed.If the system forms a truly stable emulsion (coalescence fully suppressed due to surfactant concentration above CMC in the bulk phase), C m could be replaced in the coalescence suppression formula by surface coverage corresponding to interfacial tension at CMC, as calculated from the Szysztowski equation.In other words, this would mean that the maximum practical surface coverage is not the Langmuir isotherm parameter C m , but Physical limitation or surface coverage discontinuity associated with CMC seems to be surprisingly often neglected in drop size distribution modeling when surfactants are present.Although it may be questioned how well Langmuir isotherm is performing above bulk concentrations much higher than CMC, the applied model must always be physically realistic, i.e. not predict negative surface tensions and be consistent in describing material balances in both bulk phase and at the interface.

Test cases and analysis methods
To test the proposed method, numerical simulations for some relevant test cases are carried out, and the results analyzed based on observed behavior of the drop distribution and material balances.Before presenting the actual test cases, some practical modeling aspects are discussed.
In many practical cases where surfactants are added to liquidliquid systems, the objective is to reduce drop sizes.This leads to a situation where drop sizes may vary significantly during the simulation.The present approach where surface material balance is solved along with droplet population balance may in these cases lead to a situation, where several originally significant larger drop size categories are practically useless towards the end of the simulation, i.e. there are practically no droplets of the given size, and accordingly practically no material inside the droplets or at the surface.This may lead to numerical challenges, as negligible amount of material is divided by negligible surface area.It turns out that the simulation runs much faster if those categories are excluded from the simulation after they become unused.In the present simulations, this was done by setting a threshold for the volume fraction that should remain in the largest drop size category which is still considered in the simulation.The limit was set to 0.001/NC, where NC is the original number of size classes.After reaching this low droplet volume threshold, remaining traces of drop volume and amount of material in the unused large size categories was allowed to approach zero as a first order process with a time constant set to 1000 s.Despite the fact that the amount of remaining drop volume material is very small, this residual volume and material was allocated to the largest size category still in use.This prevented oscillation between active and passive categories and kept material and volume balances exact.This approach had minimal effect to the final results, but improved time integration speed and robustness in cases where several largest size categories were unused.Although the number of active categories is reduced, due to the high order nature of the present discretization method, this is not expected to result in significant numerical error.Another option in such cases would be to adapt the whole computational grid when size categories become useless.This can be done in a straightforward manner by treating drops in the original grid as mother droplets after coalescence to be distributed into the new grid, a process for which the algorithm already exists in the present method.Practical rules or advices when the adaptation should be done is however outside the scope of this work, and thus grid adaptation is not considered here but left for the future.
One interesting further characteristic parameter not often analyzed in the population balance framework is the fraction of mixing energy used for creating new surface area, as predicted by the population balance models.This can be calculated as where rate of new surface appearing due to drop breakage (m 2 /s) is multiplied with surface tension for production rate of surface energy (W), and compared to the total mixing energy dissipated in the system.This parameter can be considered as a mixing efficiency in emulsion formation processes.All the present test cases are simulated by assuming ideally mixed batch system to simplify result analysis.V. Alopaeus Chemical Engineering Science 248 (2022) 117269

Case 1. Mixing of liquid-liquid dispersion with varying initial surfactant concentrations
In the first test case, typical physical properties for a liquid-liquid dispersion were selected.This was done to test the method performance, and to reveal interesting phenomena which may not be obvious without a model at hand, capable of predicting surfactant dynamics.Continuous phase physical properties are those of water, and dispersed phase properties are typical for organic solvents.Two surfactants are modeled, one being strongly adsorbing and the second one moderately adsorbing.These are not corresponding to any particular surfactant chemicals, but the properties are representative of the two surfactant strengths, as listed e.g. in Chang and Franses (1995).Diffusion coefficient is set to a typical value for relatively large molecules in liquid phases.Surfactant distribution coefficient between the bulk phases was set in the base case as 10 À12 , i.e. assuming a surfactant non-soluble to organic droplet phase, but additional cases were tested with a surfactant soluble only in organic phase and with equally distributing surfactant.
With this system configuration, various initial loads of surfactant were simulated.In both cases, initial surfactant loads were varied so that the surfactant concentrations ranged from far below CMC to concentrations well above it.A qualitative change in the system behavior is expected somewhere between these extremes.

Case 2. Numerical test of the proposed method with interesting surface dynamics
The second test case is mainly intended to study the numerical behavior of the proposed model.This is done by selecting one case from the previous simulations, where the system is not very dilute nor saturated.At such intermediate conditions the most interesting phenomena can often be observed, and the numerical comparisons are furthest from trivial.As there is no analytic solution to the present system, solutions with small number of categories are compared to one with a high number (100 categories), which was considered as the reference case.

Case 3. Surfactant feed in the dispersed phase and effect of phase equilibrium
In the previous cases, surfactant was assumed to be non-soluble to the dispersed organic phase.Feeding of surfactant with similar adsorption properties but which is non-soluble to aqueous phase and consequently fed with the dispersed phase was tested in this case example, as well as surfactant which is equally soluble in both phases (K = 1) but fed in the dispersed phase.In this case dynamic behavior is of primary interest, as there turned out to be interesting differences depending on the phase equilibria.Total surfactant dose is the same in each case, so concentration in the dispersed phase feed case is higher due to smaller dispersed phase volume fraction.

Case 1, medium strength surfactant
Results for Sauter mean diameter and final bulk phase concentrations ratio to CMC can be found in Fig. 5 for medium strength surfactant.Surfactant addition starts to affect Sauter mean diameter at relatively small doses, and a plateau is observed slightly above CMC doses.Final bulk liquid concentration is quite naturally increased as the initial dose increases (please note logarithmic scale in initial CMC concentration).There seems to be no reason to use much higher concentrations than CMC in this system if small droplet sizes are sought.
Fig. 6 expresses surface average final surfactant coverage.Here again it can be seen that below CMC the surface is not fully covered, and adding more surfactant increases coverage.Near the CMC, surface is close of being already maximally covered as calculated from the interfacial tension corresponding to CMC.In this case this maximum coverage is approximately 87%.With low surfactant doses, the fraction of surfactant molecules adsorbed to the actual interface is approximately 13% in this example.
Finally, energy efficiency for new surface formation as predicted by the population balance rate functions is shown in Fig. 7.The shape follows quite closely that of Sauter mean diameter.The connection, however, is indirect: as more surfactants are added, drop coalescence is first suppressed and surface tension stabilizing droplets is reduced.This leads to smaller and smaller drop sizes up to a point where they do not break anymore with significant rates, and new surface formation becomes negligible despite constant mixing power.At small surfactant doses, the system reaches steady state, but droplets continuously break and coalesce.The system is thus at steady state, but not in equilibrium, at least in the thermodynamic sense of the word.

Case 1, strong surfactant
The results for a strong surfactant are qualitatively different than for medium strength surfactant.Even in cases where initial dose is slightly above CMC, available liquid-liquid interface adsorbs most available surfactant and leaves bulk phase concentrations well below CMC if dispersed phase volume fraction is significant.It is obvious that in these cases CMC cannot guide alone proper surfactant doses, but information about the final drop surface area is needed to estimate how much surfactant is needed.
In Fig. 8, Sauter mean diameter and residual surfactant in the bulk liquid are again shown.It can be seen that in this case surfactant doses approximately 500 times CMC are needed to reach fully surfactant dominated system.For the strong surfactant, CMC is naturally much lower than for the medium strength surfactant.In terms of a true dose, around 30 times less surfactant in moles is needed in this case.However, this is not obvious from the surfactant data alone, and it depends on the final drop sizes and volume fractions.A rough estimate for excess concentration required to compensate adsorbed surfactant can be obtained by simply esti-mating amount of material needed for a single layer surface coverage predicted by the Langmuir isotherm in this case: The above equation is for continuous phase surfactant dosing.For strong surfactants and dense dispersions, this is typically higher than CMC.
Relative surface coverage for the strong surfactant is shown in Fig. 9 and fraction of mixing energy used for new surface formation in Fig. 10.These follow similar qualitative trends as for the medium surfactant case but are somewhat steeper.This indicates that changes are more rapid as surfactant is added.For surfactant doses up to 200 times of CMC, 99.9% of the added surfactant is adsorbed at the interface in this test case.

Case 2, numerical test
Although numerical tests logically should be done before analyzing the system in more details, a non-trivial test case had to  be selected first.Here the test case was selected to be the strong surfactant case where initial dose to continuous phase is 100 times larger than CMC.Total number of size categories were varied from a very small number (6 in this case) to a very high number (100), where the highest number of categories was selected as the reference case in numerical error calculation.In Fig. 11, relative errors for the final Sauter mean and bulk liquid concentrations are shown.20 categories selected earlier gives already very satisfactorily results, with relative errors being less than 0.3%.With 30 categories the numerical error is already less than 0.04%.(See Table 1) As described earlier, the number of active categories is not the same as the number of total categories in the present implementation where grid adaptation is not performed.The final number of active categories perhaps represent better the required number of variables in the simulation.Number of original and final categories are compared in Table 2.
In this case approximately 2/3 of the original categories were active at the end.It is rather surprising that the relative error is smaller than 10 % already with 6 original and 3 final active categories.These errors and number of variables correspond to those typically used in quadrature methods of moments (Alopaeus et al., 2006a).
Solution time for the selected base case (20 size categories) is under 20 s with a standard laptop computer using Matlab solver ode15s with time integration tolerances set very tight, 10 À9 .Very tight time integration tolerances were used to focus any potential errors in space discretization, which is the main objective of the present work, although much looser tolerances already produce very accurate solutions.In the tested cases, simulation time is  affected by the discontinuous nature of the system, as number of active categories changes during the simulation in a discrete manner.After each change in the number of active categories, the applied time integrator starts with shorter time steps before the time step size is again increased.Despite this, the simulation runs rapidly through several discrete changes.
In most cases slow simulation times result from challenging physical closures.In some tests with a Langmuir isotherm not limited by CMC concentration as described above, numerical challenge may result from the physical restriction of surface coverage being limited below 1.If the surface is almost fully covered, small changes in the interfacial material balances may lead to surface coverage values temporarily above 1 due to numerical reasons in time integration.In these cases, some of the applied closure models including logarithms lead to obviously non-physical complex numbers.Therefore, it is important to select the closure models so that they are not only accurate enough, but do not lead to unphysical or numerically problematic situations during normal simulations.This applies to the numerical method proposed here as well as any other similar method.This case was simulated with the strong surfactant with two separate initial doses.The first dose was 100 times CMC (one of the simulations in Case 1) when fed to the continuous phase, and 900 times CMC when fed to the dispersed phase.Since dispersed phase volume fraction was set to 0.1 in these cases, this leads to equal doses in both cases.An interesting phenomenon can be seen: Even at the steady state (end of the batch), surface concentrations depend on the drop size.This is due to the fact that the drops are still undergoing both breakage and coalescence.After each breakage process, surface area (m 2 ) increases but the total amount of surfactant (mol) in the daughter drops remain the same.This leads to lower surface concentration (mol/m 2 ) as in the original drop.The opposite happens when two droplets coalesce; the droplet born upon coalescence has a higher surface concentration than the coalescing droplets.This phenomenon is compensated by continuous mass transfer from the higher surface concentration droplets first to the continuous phase and then to the lower surface concentration droplets.When a small amount of strong surfactant is added to the system as in this test case, mass transfer is not sufficient to compensate surfactant collection to large droplets due to breakage and coalescence.In these cases, final surfactant concentrations may actually depend almost linearly on the drop sizes.When larger amounts of surfactants are added, faster mass transfer between droplet surface and continuous phase compensates this phenomenon together with suppressed droplet coalescence rates.
Surface concentrations for this test case with relatively low dose of a strong surfactant are shown in Fig. 12 for those categories with significant droplet volume (active categories).Note that despite initial surfactant doses significantly higher than CMC, the amount of surfactant is still rather small as it will be collected by the produced interface, resulting in concentrations below CMC in the bulk phases.Only with a model capable of predicting drop size dependent surface concentrations this phenomenon can be studied.
In this case, final concentrations at the surface were practically the same irrespective whether the surfactant was fed with the continuous or dispersed phase.Also, final droplet size distribution was practically similar.However, dynamical behavior was somewhat slower when the surfactant is fed to the dispersed phase.This is even more pronounced with higher surfactant doses.In Fig. 13, dynamical behavior is compared in this case for continuous phase surfactant feed and dispersed phase feed.In the case of dispersed phase feed, the system is far from steady state even after one hour of mixing.It would require several hours to reach steady state in this case.
Both feed type dynamics depend on the surfactant diffusion coefficients; with higher surfactant diffusion coefficients the   dynamical behavior is naturally also faster.However, the difference in system dynamics between the feed types seems to remain even with diffusion coefficients characteristic for smaller surfactant molecules.

Conclusions
A distributed material balance approach for population balances was extended to study surfactant dynamics in liquid-liquid systems, where material balances are formulated also for the surface phase.The method is based on high order moment conserving method of classes (HMMC) which has earlier been shown to be very accurate for drop size distributions and concentration polydispersity in the drop phase.The method was first developed in its most generic form, then relevant physical closures were discussed, and finally the method was tested with a liquid-liquid dispersion modeling with realistic physical and operation parameters.The model was fast, accurate, and was capable of capturing realistic physical behavior.Furthermore, a dispersion efficiency parameter was introduced to analyze formation of new drop surface area, based on population balance model solution.Although the present formulation is based on distributed material balance approach, reducing full 2D population balance model to pseudo-2D, there may be some special cases where full bivariate approach would be beneficial.These could be related to situations where two emulsions with similar drop size distribution but highly different concentrations are mixed.Finding these method boundaries and extending the proposed methodology to full bivariate population balances is one interesting future prospects.
The present model framework was tested for liquid-liquid dispersions, but it can be extended in a straightforward manner for gas-liquid systems and other more complex processes, such as flotation where the modeled surface variables are surfactant concentration and mass transfer dynamics, and particle attachment and detachment dynamics for various bubble sizes.It is also amenable for implementation to complex multiphase process models and computational fluid dynamics software.Due to these characteristics, it can be easily used in future closure model development, where presence of surfactants can be accounted for together with the rate function parameter optimization.

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.

Fig. 1 .
Fig. 1.Schematic representation of the model on a single droplet level.

Fig. 2 .
Fig. 2. Mass transfer and relevant equilibria in the model.

Fig. 4 .
Fig. 4. Interfacial tension as a function of surfactant equilibrium concentration with Langmuir isotherm and CMC limitation.

Fig. 5 .
Fig.5.Sauter mean diameter and final bulk concentration in medium surfactant case as functions of surfactant dose.

Fig. 6 .
Fig. 6.Final surfactant coverage as a function of initial dose for medium surfactant.

Fig. 7 .
Fig. 7. Fraction of mixing energy used for new droplet surface formation as a function of initial CMC dose.

Fig. 8 .
Fig. 8. Sauter mean diameter and final bulk concentration in strong surfactant case as functions of surfactant dose.

Fig. 9 .
Fig. 9. Final surfactant coverage as a function of initial dose for strong surfactant.

Fig. 10 .
Fig. 10.Fraction of mixing energy used for new droplet surface formation as a function of initial CMC dose for strong surfactant.

Fig. 11 .
Fig. 11.Relative errors as functions of number of categories.

Fig. 12 .
Fig. 12. Final surface concentrations as a function of drop size.

Table 1
Parameters in the emulsion test cases.
Initial drop Sauter diameter 31.4 lm 3.4.Case 3. Comparison of dispersed and continuous phase surfactant dynamics

Table 2
Comparison of number of original and final active categories.