Marangoni circulation in evaporating droplets in the presence of soluble surfactants

a r t i c l e i n f o


Introduction
The evaporation of sessile droplets is a phenomenon that has a broad practical relevance. From applications like spray cooling [1] and cleaning/drying of semiconductor surfaces [2] it can be seen that many modern technologies involve sessile droplets. In some technologies the aim is to leave a homogeneous deposition pattern of particles after evaporation of the volatile components from sessile droplets. Examples of this are inkjet printing [3][4][5], pesticide spraying [6,7] and manufacturing DNA/protein microarray slides [8].
A common issue hindering the formation of a homogeneous deposition pattern is the occurrence of the coffee-stain effect: preferential evaporation at the contact line causes a strong, outward flow resulting in a ring-like deposition pattern (given a pinned contact line) [9][10][11]. Clearly, this coffee ring is the complete opposite of a homogeneous deposition pattern. Several methods of counteracting the coffee-ring effect can be found in literature (e.g. unpinned contact lines [12], oil menisci [13], ellipsoidal particles [14]) of which an important one is what is defined here as 'Marangoni circulation'. A surface tension gradient results in an interfacial flow (the Marangoni effect), which -if the surface tension gradient is negative towards the contact line -can result in a circulating flow. This Marangoni vortex is able to suppress the coffee-stain effect [15][16][17].
The surface tension gradient required for Marangoni circulation to occur, can generally have two possible origins. First, the surface tension can become non-uniform due to thermal effects, such as evaporative cooling [18] or heated substrates [19,20]. Second, the surface tension can become non-uniform as a result of nonhomogeneous changes in composition. This happens for example during the drying of multicomponent droplets [21,22] or droplets with surfactants [23,24]. Although these causes are fairly well established, it cannot be predicted straight-forwardly whether Marangoni circulation will occur or not [25].
In this present work the focus lies on the occurrence of Marangoni circulation in evaporating droplets with soluble surfactants. It is hypothesized that using dimensionless numbers, regime plots can be composed that indicate whether Marangoni vortices can be expected or not. As a result, it may be possible to identify which parameters are relevant and what tuning they require to promote Marangoni circulation and hence a more homogeneous deposition pattern. This is done both below and above the critical micelle concentration (CMC), where surfactant monomers cluster and form micelles.
Previous numerical studies focused primarily on the evaporation of droplets without surfactants [12,21,[26][27][28] or on the evolution of nonvolatile, surfactant-laden droplets [6,[29][30][31][32][33][34][35][36][37]. Only a few numerical studies combined evaporation and surfactants. Using lubrication theory, Van Gaalen et al. [38] and Karapetsas et al. [23] both considered evaporating droplets with insoluble surfactants under partial and complete wetting conditions, respectively. Furthermore, Jung et al. [39] studied the formation of coffee rings in evaporating droplets with soluble surfactants, using a lattice gas model. However, they did not consider the underlying fluid dynamics. The present work is the first numerical study to analyze the evolution of evaporating droplets with soluble surfactants. Besides, while the number of experimental studies on evaporating droplets with surfactants is growing [8,17,[40][41][42][43][44], only limited actual flow visualisations have been made [24,45]. In this respect, the numerical work presented here is a useful addition to the small number of flow visualisations.
Results are obtained by means of a numerical model based on lubrication theory. Here, by assuming a relatively thin droplet, the Navier-Stokes equations can be simplified into a 1D height evolution equation. By combining this evolution equation with surfactant transport equations at the interface and in the bulk, velocity profiles can be calculated. Counter-intuitively to what one would expect, lubrication theory is able to capture flow topologies, such as circulation, rather well [46][47][48]. This unexpected, empirical overperformance has mathematically been validated by Krechetnikov [49].
This article is organized as follows: first, in Section 2, the numerical model is presented and explained. The lubrication equation and surfactant transport equations are introduced and the evaporation model is presented. Then, in Section 3, the numerical results are shown and analyzed. For several dimensionless parameters regime plots are composed and shown in which it is distinguished whether Marangoni circulation occurs or not. The results are compared with literature. In the final section, conclusions are drawn from the obtained results.

Mathematical model
In this section, the mathematical equations are introduced that describe the droplet evolution and surfactant transport.

Drop evolution
A droplet that is deposited on a substrate is considered. The contact angle h is assumed 'small' (h 6 40 as shown by Hu and Larson [50,51]) and the density q and dynamic viscosity l are taken constant. Also, the typical drop height H 10 À3 m, which implies that the Bond number Bo ¼ qgH 2 r ( 1 (with gravitational acceleration g and surface tension r). Thus, the effects of gravity can be disregarded [52,53]. Note that this assumption only holds, because there are no possible effects of buoyancy in a single component droplet. For multicomponent droplets, gravity cannot be neglected as shown by Edwards et al. [54] and Li et al. [55]. The Reynolds number is much smaller than unity and a cylindrical coordinate system ðr; a; zÞ is used with the assumptions of no swirl (angular velocity U a ¼ 0) and axisymmetry ( @ @a ¼ 0). Given this case, the Navier-Stokes equations can be rewritten into a 1D evolution equation for the height hðr; tÞ as a function of radial coordinate r and time t [21]: Here, p denotes the pressure inside the droplet and w e the evaporative volume flux, which is negative. From this equation it can be recognized that fluid in a higher pressure region will tend to flow towards a lower pressure region, while fluid in a lower surface tension region will tend to flow towards a higher surface tension region (the Marangoni effect). Furthermore, any evaporation of fluid will accordingly result in a decrease in local drop height. For a derivation of Eq. (1) and the velocity field (u; w), see Section S1 of the Supplementary material. The pressure in the droplet is dominated by surface tension effects. Therefore, the pressure p can be given by the Laplace pressure: Substituting the Laplace pressure in Eq. (1) shows that the droplet will tend toward a spherical cap shape. If there is no evaporation (w e ¼ 0) a droplet will tend towards an equilibrium, constant curvature shape ( @p @r ¼ @r @r ¼ 0). Note that r is also part of the derivative in Eq.
(2), because it is not necessarily constant in the presence of surfactants. Thus, rather than writing r outside the derivative as is usually done (e.g. see [12]), it should be included in it as shown by Thiele et al. [56]. Note that Eq. (2) is beyond the lubrication limit, because of the square root in the denominator. In traditional lubrication theory, this term is approximated as unity.
The boundary and initial conditions that hðr; tÞ is subjected to are given by: @h @r r¼0 ¼0; ð3Þ hðR; tÞ ¼0; ð5Þ hðr; 0Þ ¼h 0 ðrÞ: Here, R is the drop radius and h 0 is the initial drop profile, given by the previously mentioned spherical cap shape. The considered cases involve a contact line that is pinned rather than one that moves [57,58]. This means the drop radius will remain constant during evaporation, while the contact angle decreases, as opposed to a decreasing radius with a constant contact angle. Contact line pinning typically occurs for relatively small contact angles and rough substrates and is also promoted by the presence of surfactants [8,38,41,43].

Interfacial surfactant transport equation
At the liquid-air interface of the droplet a surfactant concentration Cðr; tÞ can be defined, which can be used to describe the transport of adsorbed surfactant molecules along the interface. As shown by Wong et al. [59], the surfactant transport equation is given by: Here, U t is the fluid velocity tangential to the liquid-air interface, D C the surface diffusion coefficient and J C/ is the rate with which surfactant is exchanged with the bulk through adsorption and desorption. Note that J C/ can be either positive and negative. The derivative @ @s is the surface derivative which can be written as @ @s ¼ 1 ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1þð@r hÞ 2 p @ @r . The first term on the right-hand side of Eq. (7) accounts for convective transport tangential to the interface and the second term for transport normal to the interface, which in practice boils down to a source-or sink-like effect as the interfacial curvature decreases or increases respectively. The third term denotes diffusion, the fourth term corrects for the displacement of the surface coordinates that move along the normal of the surface and the last term is the adsorption/desorption of surfactant from the bulk onto the interface and vice versa. For a derivation of Eq. (7), see Section S.2 in the Supplementary material. The surfactant concentration Cðr; tÞ is subjected to the following boundary conditions and initial condition: Cðr; 0Þ These denote the symmetry condition, no-flux condition and initial (homogeneous) surfactant concentration C 0 respectively.
Since surfactants adsorbed at the interface decrease the surface tension, an equation of state for r is required to close the problem.
In this work, the Frumkin equation is used, which is analogous to the Langmuir isotherm and considers that surfactant molecules occupy a finite amount of space at the interface [60]. The Frumkin equation is given by: Here, r 0 is the surface tension of the pure liquid, R u is the universal gas constant, T the temperature (which is assumed constant) and

Bulk surfactant transport equation
Besides surfactant being adsorbed onto the interface, there is also surfactant dissolved in the bulk of the droplet. The bulk monomer concentration / is considered to be a function of r and t and independent of z, which is allowed as rapid vertical diffusion is assumed (as outlined by [35,61]). In order to let the concentration be independent of hðr; tÞ, the evolution of the bulk monomer concentration is described as a function of wðr; tÞ ¼ /ðr; tÞhðr; tÞ, as introduced by Thiele et al. [62]. The transport equation can then be written as: In this equation, three terms can be distinguished. The first term is a transport term (with D / the diffusion coefficient of surfactant monomers in the bulk), which takes into account convective and diffusive transport, the second term is the adsorption/desorption of surfactant from the bulk onto the interface and vice versa, including a factor that compensates for the interface geometry, and the third term is the micelle formation rate J M/ multiplied with the preferred number of surfactant monomers N which form a single micelle.
When the bulk surfactant concentration exceeds a certain threshold, it becomes energetically more favorable for the molecules to cluster together and form micelles. This threshold is called the critical micelle concentration (CMC). In practice, the number of surfactant monomers that form a single micelle tends to strongly prefer a single value N, which is also assumed in this work [63].
Similarly to the monomer bulk concentration /, the micelle bulk concentration M is given in terms of fðr; tÞ ¼ Mðr; tÞhðr; tÞ: Here, D M is the diffusion coefficient of micelles. Note that it is assumed here that micelles do not adsorb directly onto the interface. This is a common assumption in literature [30,33,35], that follows from the fact that surfactant monomers need to dissociate from micelles before they can be adsorbed onto the interface [64].
Since the micelle formation/decomposition process is modelled as a single step, it is consistent to only allow individual monomers to adsorb onto the interface. Future models, however, could include multi-step models. The bulk concentrations wðr; tÞ and fðr; tÞ are subject to the following boundary conditions and initial conditions: wðr; 0Þ ¼/ 0 hðr; tÞ; ð16Þ fðr; 0Þ ¼M 0 hðr; tÞ: Similarly to the boundary and initial conditions of C, these denote the symmetry condition, the no-flux condition at the contact line and the initial, constant bulk concentrations / 0 and M 0 respectively. All initial surfactant concentrations (C 0 ; / 0 and M 0 ) are always chosen to be in equilibrium, so initially J C/ ¼ J M/ ¼ 0. A derivation of both Eqs. (12) and (13) is given in Section S.3 of the Supplementary material.

Surfactant adsorption
The transport between the interface and the bulk is a continuous adsorption and desorption of molecules that overall tends towards a dynamic equilibrium. Furthermore, there is a limited amount of space available at the interface for surfactants, so as the interfacial concentration increases the adsorption rate tends to decrease, while the desorption rate increases. This behavior can be described by the following reaction equation [33,35]: Here, S (¼ 1 À C C1 ) indicates the fraction of available space at the interface, k C a is the interfacial adsorption coefficient, and k C d the interfacial desorption coefficient. Thus, the interfacial adsorption flux is given by: In this equation it can indeed be recognized that the first, adsorption term increases as the bulk concentration increases and approaches zero as C ! C 1 , while the second, desorption term becomes more negative as C increases. This behavior corresponds indeed to Eq. (18). Note that surfactant adsorption onto the substrate is not taken under consideration here, because this will not have a direct effect on the internal flow. Marangoni flow can only occur as a result of a surface tension gradient at a free interface. Any indirect influences of substrate sorption on the flow -surfactant being removed or added to the bulk -are considered less significant. Similar to Eq. (18), a reaction equation can be written for the formation of micelles: Here, k M a is the micelle formation coefficient and k M d the micelle decomposition coefficient. This reaction equation can subsequently be written into a micelle formation rate J M/ as well [30,33,35]: Given this equation, it is possible to approximate the reaction constants in terms of the CMC and concentrations. It can be seen that the equilibrium concentration is given by Uð¼ NM þ /Þ is the total concentration of surfactant monomers in the bulk, substitution of the initial total concentration U 0 results in: Here, it should be noted that in equilibrium / ¼ CMC, given that there are micelles present. Of course, it is still required to choose one of the reaction constants to define the time scales of the reactions.

Evaporation
The evaporative volume flux w e is given by: Here, D v is the vapor diffusivity in air, M l the molar mass of the liquid, p sat;l the liquid saturation pressure, which is assumed constant, , with p l the local vapor pressure. The derivative @ @n is the spatial derivative along the normal vector (pointing outwards) [21].
In order to calculate @p l @n , it is required to describe the vapor field. As shown by Deegan [11] and Hu and Larson [65], vapor diffusion can be considered instantaneous. Together with the assumption of no convection, this implies that the vapor field can be described by the Laplace equation: The corresponding boundary conditions are given by: Here RH is the relative humidity. These equations represent saturated vapor at the droplet surface, no vapor penetration at the substrate and ambient relative humidity far away from the droplet. An analytical solution to Eq. (24) is derived in Section S.4 of the Supplementary material.

Dimensionless parameters
In order to analyze when Marangoni circulation occurs, it is helpful to define dimensionless parameters that describe the droplet characteristics. To this purpose, Table 1 can be compiled, which gives the ratio between several important scales. Typical numerical values of variables used for the dimensionless numbers are given in Section S.6 of the Supplementary material.
One may note that the traditional capillarity number Ca ¼ lw e =r 0 is not listed in Table 1. This dimensionless number is omitted, because for any physically relevant case Ca ( 1. Varying Ca therefore has negligible effect on the drop evolution and internal flow field, because surface tension forces always dominate over viscous drag forces, resulting in little deviations from the spherical cap shape. Simulations show that changing Ca by a factor 10 3 only results in a 0.0067% change in total pressure gradient from the drop center toward the contact line.

Results and discussion
In this section, the effect of changes in several dimensionless parameters on the occurrence of Marangoni circulation during droplet evaporation with soluble surfactants is analyzed. First, cases without micelles and diffusion are analyzed and subsequently surface and bulk diffusion are taken into account. After that, cases with micelles are considered.

Below CMC without diffusion
Simulations are carried out for various ranges of the dimensionless parameters. The numerical procedure is outlined in Section S.5 of the Supplementary material. The initial contact angle h 0 is always set to 25 and the contact line is pinned. During each simulation, at least 200 time steps are calculated ($5% of the drying time), after which the velocity field is analyzed.
Only the initial stage of the evaporation process is considered, because this allows one to capture a 'snapshot' of the flow dynamics. It is still required to allow the internal flow dynamics to evolve sufficiently, ensuring that the initial state has no significant effect on the flow field. If the full evaporation process would be taken into account, this would give a skewed perspective on the flow dynamics, because the considered dimensionless numbers change over time. For example, the height decreases over time and the interfacial surfactant concentration increases over time. Also, the bulk concentration may exceed the CMC at some point in time. Analysis of the entire evaporation process is outside the scope of this paper and can be considered in future work. Nevertheless, the results of this work can still be used to predict flow patterns during the entire evaporation process, as long as the relevant dimensionless numbers can be estimated.
Visual analysis of the velocity field has led to the definition of three separate regimes: the 'coffee-ring regime', where the flow field looks similar to pure droplet evaporation, the 'Marangoni regime', where a clear Marangoni eddy can be distinguished, and the 'transition regime', which shows behavior of both the coffeering and Marangoni regime. Representative velocity plots of the three regimes are shown in Fig. 1. In all three velocity regimes an outward, capillary flow can be distinguished, that is caused by selective evaporation at the contact line. However, for the Marangoni regime (and in some degree for the transition regime) there is a backflow close to the interface that convects adsorbed surfactant towards the center, where it desorbs into the bulk again. All three regimes have a zero-fluid velocity at both the liquid-air and liquid-solid interface at r = 0, which is in accordance with the 'Hairy ball theorem' (Poincaré-Brouwer theorem). This theorem predicts that there necessarily exists at least one zero-velocity point on the surfaces of compressible liquids and interfaces allowing mass, energy and momentum transport [66].
Calculating the mean, dimensionless vorticity x over the middle area where the Marangoni vortex tends to appear (around R=2 < r < 5=6R) shows that x 6 À1:2 Á 10 À4 corresponds to the Marangoni regime, À1:2 Á 10 À4 < x 6 À0:8 Á 10 À4 to the transition regime and x > À0:8 Á 10 À4 to the coffee-ring regime. Here, the time scale t d ¼ R 2 =D v is used to make the vorticity dimensionless.
Alternatively, as a cross-validation the resulting velocity fields can be classified by calculating the fraction of velocity vectors opposing the typical coffee-ring flow. Here, numerical analysis shows that the Marangoni regime corresponds to more than 24% of the radial velocity vectors and more than 9% of the axial velocity vectors pointing to the center and upward respectively. Similar, the transition regime corresponds to velocity fields that are not in the Marangoni regime with more than 22% of the radial velocity vectors and more than 6% of the axial velocity vectors pointing to the center and upward respectively. If a lower proportion of the velocity vectors points centerwise or upwards, the velocity field is considered to be in the coffee-ring regime. This gives similar results to the classification with vorticity and small variations of the values do not yield significantly different classification results. Further classification of regimes in this paper is made through the mean vorticity with frequent, random visual checks for extra verification. Fig. 2 shows the regime charts for simulations without diffusion and without micelles. As can be seen in Fig. 2a, for Tr < 1, both 1=Ev and Tr have very similar effects on the occurence of Marangoni flow. This makes sense, because as 1=Ev decreases, the evaporative velocity starts to dominate over the Marangoni velocity. The evaporative effects are thus too strong for the Marangoni effect to counter, resulting in a coffee-ring regime flow. Similarly, as Tr decreases, the effects of adsorption and desorption become less significant (given that De remains constant). This will result in behavior as if the surfactant is insoluble, which is described by the coffee-ring regime [38]. Because the relative strength of the Marangoni effect increases, the interfacial concentration is kept homogeneous since the fluid velocity close to the interface is reduced and the concentration increase at the drop apex due to curvature effects becomes more dominant. At the same time less adsorption is taking place, so the increased bulk concentration close to the contact line will not result in enough interfacial adsorption to cause a positive concentration gradient. That increasing Tr tends to promote Marangoni circulation is also found by Jung et al., who used a lattice-gas model [39].
For Tr > 1, the occurrence of Marangoni circulation becomes mostly limited by 1=Ev. Only if a certain threshold value of 1=Ev is exceeded, Marangoni vortices can be formed. At the same time, Tr is largely irrelevant, because the desorption kinetics are no longer dominated by the Marangoni velocity and thus are not a limiting factor anymore. That 1=Ev still is a limiting factor however, is not surprising. Decreasing surfactant concentration C 0 or increasing the evaporation rate (e.g. through an increase in D v ) will eventually always result in a coffee-ring flow, because either the droplet can be considered pure or the evaporation becomes too dominant. It is important to note that a correct interpretation of Tr involves not only the desorption rate, but also the adsorption rate. If Tr is varied by changing k d the value of k a changes accordingly, given that De should remain constant. Increasing Tr thus does not only mean that interfacial desorption of surfactants becomes more dominant with respect to the Marangoni velocity, but interfacial sorption in general becomes more dominant.
Another relevant observation is that the dimensionless numbers do not necessarily give information about the strength of the Marangoni vortex. As an illustration, Tr can be modified by changing C 0 and by changing k d . An increase in C 0 results in a proportional increase in the absolute velocity however, while a decrease in k d will generally not have that influence. This becomes especially relevant when other physical effects are involved, such as thermal Marangoni flow or the deposition of a solute. For example, Jung et al. conclude that increasing the initial surfactant concentration tends to increase the Marangoni strength, thus suppressing the formation of coffee-ring deposits [39]. With respect to Fig. 2a this can only be explained by also considering the absolute velocity rather than only relative to other time scales.
In Fig. 2b the influence of the Desorption number De is shown. An increase of Tr still tends to promote Marangoni circulation, but this only holds approximately for De > 10. As the value of De decreases, the contribution of adsorption becomes increasingly similar in magnitude to desorption. As a result, Marangoni circulation is increasingly suppressed. The reason for this suppression is that De effectively is the solubility of the surfactant. As De decreases, the surfactant solubility also decreases and the surfactant increasingly starts to behave as insoluble. This effectively means that the bulk concentration becomes relatively small with respect to the interfacial concentration. As a consequence, adsorption onto the interface, resulting from concentration increases in the bulk, will become less significant. Therefore, positive interfacial concentration gradients towards the contact line will no longer arise. On the contrary, the interfacial concentration gradient will tend to be negative towards the contact line as a result of the interface shrinking fastest at the drop apex. The result is a flow in the coffee-ring regime. The flow close to the interface is here reduced by the Marangoni effect to counter any positive surfactant gradients. This flow behavior was previously found by Van Gaalen et al. [38] for insoluble surfactants. The behavior of Fig. 2a and 2b is reflected in Fig. 2c: if evaporation becomes more dominant the flow tends to transition to the coffee-ring regime. Similar, as the desorption number is decreased Marangoni circulation will not appear anymore some point.
The results in Fig. 2 are consistent with experiments performed by Marin et al. [24]. They performed experiments both below and above CMC with the surfactants polysorbate 80 (P80) and sodium dodecyl sulfate (SDS). P80, is a large, slow surfactant, with low solubility (CMC = 0.012 mM [67]), while SDS is a small, fast surfactant, with higher solubility (CMC = 8.2 mM [68]).
For P80 Marin et al. reported a suppression of thermal Marangoni flow and a severe reduction of the interfacial flow strength. This is also predicted by the numerical model, given the physical properties of P80. Slow adsorption kinetics correspond to a low Tr and low solubility to a low De. In all three subfigures of Fig. 2 it can be seen that these low values would predict a coffee-ring flow and not a solutal Marangoni flow. Furthermore, a low De also results in a reduced interfacial velocity and thus suppression of any thermal Marangoni flow, which is consistent with the experimental results.
For SDS, on the other hand, Marin et al. reported completely opposite behavior. SDS tends to increase the strength of Marangoni circulation. The fast adsorption kinetics of SDS correspond to a high Tr and the high solubility to a high De. This would place SDS in the Marangoni regime, as can be seen in Fig. 2. Consequently, it is not surprising to see an increase in Marangoni circulation, since SDS increases the magnitude of the already negative thermal surface tension gradient.
From these two examples it becomes clear how two different surfactants can have opposite effects on the flow in drops and thus it underlines the explanatory power of the numerical model.

Below CMC with diffusion
Given the results in SubSection 3.1, additional degrees of freedom can be introduced by setting the bulk and surface diffusion coefficients to a nonzero value. In this way, two different regime plots can be drawn: one with surface diffusion, but without bulk diffusion, and one with bulk diffusion, but without surface diffusion. The result are displayed in Fig. 3.
As can be seen in Fig. 3a, decreasing the surface Péclet number tends to suppress Marangoni circulation. This makes sense, because as the effect of surface diffusion increases, a smaller Marangoni flow is required to counter any interfacial concentration gradient. Adsorption from the bulk causes a surface tension gradient, but this gradient is increasingly counteracted by surfactant diffusion at the interface. The result tends increasingly towards a coffee-ring-regime flow as Pe s decreases to and becomes smaller than unity.
Similar behavior can be distinguished when the bulk Péclet number is varied. As Pe b is decreased, Marangoni circulation tends to be suppressed, because the bulk concentration increases less sharply. This way, adsorption of surfactant from the bulk onto the interface occurs over a larger area at the contact line, resulting in a flattening of the surface tension gradient. Since at the same time the surfactant concentration increases towards the drop apex due to the shrinking of available interface, a nearly constant interfacial surfactant concentration arises. Thus, the flow will tend towards the coffee-ring regime as the influence of bulk diffusion increases. Interestingly enough, the transition from the Marangoni regime to the coffee-ring regime already starts to occur for Pe b < 10 5 , which is a factor 10 4 higher than the surface Péclet number at which this transition happens. This implies that bulk diffusion has a larger influence than surface diffusion.

Above CMC without monomer diffusion
As the bulk concentration increases, at some point it becomes energetically more favorable for surfactant molecules to cluster together in the form of micelles rather than as separate monomers. This will influence the internal droplet dynamics, since as the local monomer concentration increases in the bulk, only part of the surfactant monomers will adsorb onto the interface and another part will form micelles. In Fig. 4 the effect of micelles on the internal flow patterns is shown. In the simulations the initial bulk concentration starts at CMC. Fig. 4a shows that as the micelle transport number Tr M increases, and thus the rates with which micelles are formed and decompose increase, Marangoni circulation seems to become increasingly suppressed. This suppression can be explained by the fact that the formation of micelles reduces the absorption of bulk monomers onto the interface, because the bulk concentration is now also reduced through an additional mechanism. Both adsorption onto the interface and formation of micelles can now reduce a high bulk concentration. Since less surfactant is absorbed onto the interface than without micelle formation, the surface tension gradient becomes less pronounced, which increasingly results in suppression of Marangoni circulation.
Regarding the micelle decomposition number, it is shown in Fig. 4b that for De M < 1 the value of De M is largely irrelevant.
Here, it is used that c 0 ¼ CMC. Given this reformulation it becomes clear that varying the initial micelle concentration M 0 has a proportional effect on De M and vice versa. A high value of De M (e.g. larger than unity) would thus correspond to a relatively low value of M 0 , which implies that the shift in Fig. 4b is a transitional effect: the flow dynamics have not fully transitioned from below CMC to above CMC for De M > 1. In Fig. 4c it can be seen that the micelle Péclet number Pe M has a similar effect on the flow dynamics as Pe s and Pe b . As Pe M is reduced, the transition from coffee-ring regime to Marangoni regime shifts to a higher Tr. This is caused by micelles being transported inward as a result of diffusion. There, they decompose into surfactant monomers, effectively reducing the bulk concentration gradient. Subsequently, the interfacial concentration gradient is also reduced, counteracting Marangoni circulation.
To summarize, the simulations predict that for concentrations above CMC Marangoni circulation becomes more suppressed than below CMC, although this effect is minor. In experiments however, the influence of surfactants tends to increase even more beyond the CMC. For example, Marin et al. [24] reported that for experiments with P80 above CMC the surface velocity is even more reduced than below CMC and for very high concentrations even reversed. Furthermore, they show that for SDS the Marangoni circulation becomes even stronger above CMC than below CMC and even report several Marangoni vortices at once. Similarly, Sempels et al. [45] reported for Triton X-100 an increased strength of the Marangoni flow as the surfactant concentration is increased, while it is already above the CMC.
These differences between simulations and experiments can possibly be explained by physical effects that have not been taken into account in the numerical model. For example, micelles may adsorb directly onto the interface, without the need to first decompose into bulk monomers. Models used in literature (including this work) usually assume a single step monomer adsorption model [30,33,35], but in reality the adsorption process is more complex, especially when micelles are involved [64]. This is even more the case for high deviations from equilibrium. Direct adsorption of micelles would indeed explain why the experiments show an increased Marangoni flow above CMC: for fast surfactants (e.g. SDS) there is more adsorption/desorption resulting in higher surface tension gradients.
However, this would not explain the even further reduced and reversed surface velocity for P80 and the dual vortices for SDS. This may possibly be attributed to other adsorption and transport effects, such as adsorption onto the substrate (or even micelle formation on the substrate [69]) and transport of surfactant from the substrate to the interface through the contact line [35]. This may have significant influence on the flow.
An alternative factor that could play a role in the experiments, is the influence of micelles on the fluid properties of the droplet. As an illustration, it is well known that particles (like micelles) tend to increase the viscosity of a fluid [70,71]. Especially if larger, more complex aggregates are formed [72], this may play a significant role. Also, the shape of the particles can play a major role. For example, while spherical particles do not tend to influence the flow significantly, ellipsoidal particles tend to aggregate at the interface in loosely packed structures [73][74][75][76]. This can either increase [14] or decrease [77] the surface tension and keep the particles from flowing towards the contact line. A possible way of modeling this, would be to allow adsorption-like behavior for the micelles, combined with an equation of state for surface tension. The self-assembly could then be modeled through a concentration-dependent resistance to flow and increased 'adsorption' with surface concentration.
Nevertheless, it can be concluded that while the model is consistent with experiments below CMC, it deviates above CMC. This shows the need for models that capture surfactant kinetics above CMC in a more detailed way than the standard models that are used in this work and in literature [30,33,35].

Conclusion
A numerical model was developed to predict the flow in evaporating droplets with soluble surfactants. The drop evolution and flow behavior have been modeled with lubrication theory and the surfactants were implemented by means of coupled convection-diffusionadsorption equations. Simulations were carried out for variations in several dimensionless parameters, over a broad range. The discovered changes in flow characteristics were compared with both experimental and numerical results from literature.
Below CMC, three parameters (Tr; Ev; De) were analyzed for cases without diffusion and two additional parameters (Pe s ; Pe b ) were considered for cases with diffusion. The effects of these parameters on the internal flow patterns were explained and found to be consistent with experimental and numerical results from literature [24,39].
Above CMC, three additional parameters (Tr M ; De M ; Pe M ) were analyzed. Although these results could be explained intuitively, they were found to differ from experimental results from literature. In experiments, the effect of surfactants on the flow properties tends to become increasingly more dominant as the surfactant concentration increases and even results in different qualitative behavior [24,45], while in the simulations the influence of micelles tends to be rather small. These differences can possibly be explained by more complex adsorption and transport mechanisms, such as micelle formation on the substrate [69], and fluid properties changing due to micelles, both effects that have not been taken into account in the model. This shows that more detailed models are needed to capture all relevant dynamics of micelles. These effects are not considered in current state-of-the-art models [30,33,35].
Nevertheless, the results agree with the hypothesis that was made, namely that using dimensionless numbers regime plots can be drawn to predict whether Marangoni vortices will arise. This was done both below and above CMC and the qualitative agreement with experiments is quite good below CMC.
The relevance of these findings lies first of all in the ability to explain experimental results. The numerical model shows why surfactants can have opposite influences on flow dynamics in evaporating droplets. Furthermore, this paper shows the predictive power of the model. The results can be used to predict and understand the flow dynamics in a surfactant-laden droplet. For example, surfactants that are larger and slower than P80 can be assumed to reduce the interfacial flow, while surfactants that are smaller and faster than SDS will only accelerate the Marangoni circulation. Also, it is shown that increasing the surfactant concentration will generally not increase the likelihood of encountering flow circulation (except for very low concentrations as can be seen in Fig. 2a). This counter-intuitive phenomenon is also implied by the results of Marin et al., where they show that an increase in P80 only slows the interfacial velocity further.
The explanatory and predictive power of the model is relevant in a broad range of applications. Evaporating sessile droplets are applied in various technologies, such as inkjet printing [3,4], pesticides [6,7], surface patterning [78] and spray cooling [1]. Close to all of these technologies involve surfactants, either on purpose or through contamination. It is thus of crucial importance to have a deep understanding of the effect of this component, in order to control the internal flow and deposition pattern.
Future research opportunities lie in expansion of the micelle model. The current model [30,33,35] is not able to fully explain the experimental results [24,45] and thus additional physical effects need to be added. For example, direct adsorption of micelles onto the interface and monomer adsorption onto the substrate may be able to explain experimental results above CMC. Furthermore, it may also be useful to investigate the full evaporation process of the droplet, including moving contact lines, to be able to analyze drying patterns. This will enable the control of the final deposition in technologies as inkjet printing.

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.