## Journal of Colloid and Interface Science

胶体与界面科学杂志

第584卷，2021年2月15日，第622-633页

# Regular Article 常规文章Marangoni circulation in evaporating droplets in the presence of soluble surfactants

在可溶性表面活性剂存在下蒸发液滴中的马兰戈尼循环

## Graphical abstract 图形摘要

## Keywords 关键字

液滴可溶性表面活性剂蒸发润滑近似马兰戈尼流量

## 1. Introduction 1. 引言

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].

无柄液滴的蒸发是一种具有广泛实际意义的现象。从喷雾冷却[1]和半导体表面的清洁/干燥[2]等应用中可以看出，许多现代技术都涉及无柄液滴。在一些技术中，目的是在无柄液滴中的挥发性成分蒸发后留下均匀的颗粒沉积模式。这方面的例子包括喷墨打印[3]、[4]、[5]、杀虫剂喷洒[6]、[7]和制造DNA/蛋白质微阵列载玻片[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].

阻碍均匀沉积模式形成的一个常见问题是咖啡渍效应的发生：接触线处的优先蒸发导致强烈的向外流动，从而产生环状沉积模式（给定固定接触线）[9]、[10]、[11]。显然，这个咖啡环与均匀的沉积模式完全相反。在文献中可以找到几种抵消咖啡环效应的方法（例如，无钉接触线[12]、油半月板[13]、椭球体颗粒[14]），其中重要的一种是这里定义的“马兰戈尼环流”。表面张力梯度导致界面流动（马兰戈尼效应），如果表面张力梯度朝接触线方向为负，则会导致循环流动。这种马兰戈尼涡旋能够抑制咖啡渍效应[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 non-homogeneous 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].

马兰戈尼环流发生所需的表面张力梯度通常有两个可能的起源。首先，由于热效应，如蒸发冷却[18]或加热基材[19]，[20]，表面张力可能变得不均匀。其次，由于成分的非均匀变化，表面张力会变得不均匀。例如，在多组分液滴[21]、[22]或含有表面活性剂[23]、[24]的液滴干燥过程中，就会发生这种情况。尽管这些病因已经相当明确，但无法直接预测是否会发生马兰戈尼环流[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.

在本研究中，重点在于用可溶性表面活性剂蒸发液滴中马兰戈尼循环的发生。据推测，使用无量纲数，可以组成表明是否可以预期的马兰戈尼涡旋的政权图。因此，有可能确定哪些参数是相关的，以及它们需要什么调整来促进马兰戈尼环流，从而形成更均匀的沉积模式。这是在临界胶束浓度 （CMC） 以下和上方完成的，其中表面活性剂单体聚集并形成胶束。

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.

以往的数值研究主要集中在不含表面活性剂的液滴的蒸发[12]、[21]、[26]、[27]、[28]或非挥发性、富含表面活性剂的液滴的演变[6]、[29]、[30]、[31]、[32]、[33]、[34]、[35]、[36]、[37]。只有少数数值研究将蒸发和表面活性剂结合起来。Van Gaalen等[38]和Karapetsas等[23]分别利用润滑理论考虑了在部分润湿和完全润湿条件下蒸发含有不溶性表面活性剂的液滴。此外，Jung等[39]使用晶格气体模型研究了含有可溶性表面活性剂的蒸发液滴中咖啡环的形成。然而，他们没有考虑潜在的流体动力学。本研究是首次用可溶性表面活性剂分析蒸发液滴演变的数值研究。此外，虽然关于用表面活性剂蒸发液滴的实验研究数量正在增加[8]、[17]、[40]、[41]、[42]、[43]、[44]，但只有有限的实际流动可视化[24]，[45]。在这方面，这里介绍的数值工作是对少量流动可视化的有用补充。

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.

## 2. Mathematical model

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

### 2.1. Drop evolution

A droplet that is deposited on a substrate is considered. The contact angle $\theta $ is assumed ‘small’ ($\theta \u2a7d{40}^{\circ}$ as shown by Hu and Larson [50], [51]) and the density $\rho $ and dynamic viscosity $\mu $ are taken constant. Also, the typical drop height $H\le {10}^{-3}$ m, which implies that the Bond number $\mathit{Bo}=\frac{\rho {\mathit{gH}}^{2}}{\sigma}\ll 1$ (with gravitational acceleration *g* and surface tension $\sigma $). 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,\alpha ,z)$ is used with the assumptions of no swirl (angular velocity ${U}_{\alpha}=0$) and axisymmetry ($\frac{\partial}{\partial \alpha}=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]:(1)$\frac{\partial h}{\partial t}=\frac{1}{r\mu}\frac{\partial}{\partial r}\left(\frac{{\mathit{rh}}^{3}}{3}\frac{\partial p}{\partial r}-\frac{{\mathit{rh}}^{2}}{2}\frac{\partial \sigma}{\partial r}\right)+{w}_{e}\text{.}$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:(2)$p=-\frac{1}{r}\frac{\partial}{\partial r}\left(\sigma \frac{r{\partial}_{r}h}{\sqrt{1+{({\partial}_{r}h)}^{2}}}\right)\text{.}$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 ($\frac{\partial p}{\partial r}=\frac{\partial \sigma}{\partial r}=0$). Note that $\sigma $ is also part of the derivative in Eq. (2), because it is not necessarily constant in the presence of surfactants. Thus, rather than writing $\sigma $ 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:(3)${\left(\frac{\partial h}{\partial r}\right)}_{r=0}=0,$(4)${\left(\frac{{\partial}^{3}h}{\partial {r}^{3}}\right)}_{r=0}=0,$(5)$h(R,t)=0,$(6)$h(r,0)={h}_{0}(r)\text{.}$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].

### 2.2. Interfacial surfactant transport equation

At the liquid–air interface of the droplet a surfactant concentration $\mathrm{\Gamma}(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:(7)$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Gamma}}{\partial t}=-& \frac{1}{r}\frac{\partial ({\mathit{rU}}_{t}\mathrm{\Gamma})}{\partial s}+\frac{\mathrm{\Gamma}{\partial}_{t}h}{1+{({\partial}_{r}h)}^{2}}\left(\frac{{({\partial}_{r}h)}^{2}}{1+{({\partial}_{r}h)}^{2}}+\frac{1}{r}\frac{\partial h}{\partial r}\right)+\hfill \\ \hfill & \frac{{D}_{\mathrm{\Gamma}}}{r}\frac{\partial (r{\partial}_{s}\mathrm{\Gamma})}{\partial s}+\frac{\partial h}{\partial t}\frac{{\partial}_{r}h}{1+{({\partial}_{r}h)}^{2}}\frac{\partial \mathrm{\Gamma}}{\partial r}+{J}_{\mathrm{\Gamma}\varphi}\text{.}\hfill \end{array}$Here, ${U}_{t}$ is the fluid velocity tangential to the liquid–air interface, ${D}_{\mathrm{\Gamma}}$ the surface diffusion coefficient and ${J}_{\mathrm{\Gamma}\varphi}$ is the rate with which surfactant is exchanged with the bulk through adsorption and desorption. Note that ${J}_{\mathrm{\Gamma}\varphi}$ can be either positive and negative. The derivative $\frac{\partial}{\partial s}$ is the surface derivative which can be written as $\frac{\partial}{\partial s}=\frac{1}{\sqrt{1+{({\partial}_{r}h)}^{2}}}\frac{\partial}{\partial 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 $\mathrm{\Gamma}(r,t)$ is subjected to the following boundary conditions and initial condition:(8)${\left(\frac{\partial \mathrm{\Gamma}}{\partial r}\right)}_{r=0}=0,$(9)${\left(\frac{\partial \mathrm{\Gamma}}{\partial r}\right)}_{r=R}=0,$(10)$\mathrm{\Gamma}(r,0)={\mathrm{\Gamma}}_{0}\text{.}$These denote the symmetry condition, no-flux condition and initial (homogeneous) surfactant concentration ${\mathrm{\Gamma}}_{0}$ respectively.

Since surfactants adsorbed at the interface decrease the surface tension, an equation of state for $\sigma $ 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:(11)$\sigma ={\sigma}_{0}+{R}_{u}T{\mathrm{\Gamma}}_{\infty}\mathrm{ln}\left(1-\frac{\mathrm{\Gamma}}{{\mathrm{\Gamma}}_{\infty}}\right)\text{.}$Here, ${\sigma}_{0}$ is the surface tension of the pure liquid, ${R}_{u}$ is the universal gas constant, *T* the temperature (which is assumed constant) and ${\mathrm{\Gamma}}_{\infty}$ the maximum possible surfactant concentration. Note that for $\mathrm{\Gamma}\ll {\mathrm{\Gamma}}_{\infty}$ the equation can be approximated by a linear, dilute equation of state.

Given this equation of state for surface tension and the radial velocity (as defined by Equation (S.7) in the Supplementary material) a typical Marangoni velocity can be defined as ${u}_{\mathit{Mar}}=\frac{{\mathit{HR}}_{u}T{\mathrm{\Gamma}}_{0}}{\mu R}$. This typical velocity is used in SubSection 2.6 to define several relevant dimensionless numbers.

### 2.3. 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 $\varphi $ 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 $\psi (r,t)=\varphi (r,t)h(r,t)$, as introduced by Thiele et al. [62]. The transport equation can then be written as:(12)$\frac{\partial \psi}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\left(\frac{{\mathit{rh}}^{2}\psi}{3\mu}\frac{\partial p}{\partial r}-\frac{\mathit{rh}\psi}{2\mu}\frac{\partial \sigma}{\partial r}+{D}_{\varphi}h\frac{\partial \varphi}{\partial r}\right)-{J}_{\mathrm{\Gamma}\varphi}\sqrt{1+{({\partial}_{r}h)}^{2}}-{J}_{M\varphi}N\text{.}$In this equation, three terms can be distinguished. The first term is a transport term (with ${D}_{\varphi}$ 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\varphi}$ 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 $\varphi $, the micelle bulk concentration *M* is given in terms of $\zeta (r,t)=M(r,t)h(r,t)$:(13)$\frac{\partial \zeta}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\left(\frac{{\mathit{rh}}^{2}\zeta}{3\mu}\frac{\partial p}{\partial r}-\frac{\mathit{rh}\zeta}{2\mu}\frac{\partial \sigma}{\partial r}+{D}_{M}h\frac{\partial M}{\partial r}\right)+{J}_{M\varphi}\text{.}$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 $\psi (r,t)$ and $\zeta (r,t)$ are subject to the following boundary conditions and initial conditions:(14)${\left(\frac{\partial \psi}{\partial r}\right)}_{r=0}={\left(\frac{\partial \zeta}{\partial r}\right)}_{r=0}=0,$(15)${\left(\frac{\partial \psi}{\partial r}\right)}_{r=R}={\left(\frac{\partial \zeta}{\partial r}\right)}_{r=R}=0,$(16)$\psi (r,0)={\varphi}_{0}h(r,t),$(17)$\zeta (r,0)={M}_{0}h(r,t)\text{.}$Similarly to the boundary and initial conditions of $\mathrm{\Gamma}$, these denote the symmetry condition, the no-flux condition at the contact line and the initial, constant bulk concentrations ${\varphi}_{0}$ and ${M}_{0}$ respectively. All initial surfactant concentrations (${\mathrm{\Gamma}}_{0},{\varphi}_{0}$ and ${M}_{0}$) are always chosen to be in equilibrium, so initially ${J}_{\mathrm{\Gamma}\varphi}={J}_{M\varphi}=0$. A derivation of both Eqs. (12), (13) is given in Section S.3 of the Supplementary material.

### 2.4. 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]:(18)$S+\varphi \underset{\phantom{\rule{0.12em}{0ex}}{k}_{a}^{\mathrm{\Gamma}}\phantom{\rule{0.12em}{0ex}}}{\overset{{k}_{d}^{\mathrm{\Gamma}}}{\rightleftharpoons}}\mathrm{\Gamma}\text{.}$Here, *S* ($=1-\frac{\mathrm{\Gamma}}{{\mathrm{\Gamma}}_{\infty}}$) indicates the fraction of available space at the interface, ${k}_{a}^{\mathrm{\Gamma}}$ is the interfacial adsorption coefficient, and ${k}_{d}^{\mathrm{\Gamma}}$ the interfacial desorption coefficient. Thus, the interfacial adsorption flux is given by:(19)${J}_{\mathrm{\Gamma}\varphi}={k}_{a}^{\mathrm{\Gamma}}\varphi \left(1-\frac{\mathrm{\Gamma}}{{\mathrm{\Gamma}}_{\infty}}\right)-{k}_{d}^{\mathrm{\Gamma}}\mathrm{\Gamma}\text{.}$In this equation it can indeed be recognized that the first, adsorption term increases as the bulk concentration increases and approaches zero as $\mathrm{\Gamma}\to {\mathrm{\Gamma}}_{\infty}$, while the second, desorption term becomes more negative as $\mathrm{\Gamma}$ 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:(20)$N\varphi \underset{\phantom{\rule{0.12em}{0ex}}{k}_{a}^{M}\phantom{\rule{0.12em}{0ex}}}{\overset{{k}_{d}^{M}}{\rightleftharpoons}}M\text{.}$Here, ${k}_{a}^{M}$ is the micelle formation coefficient and ${k}_{d}^{M}$ the micelle decomposition coefficient. This reaction equation can subsequently be written into a micelle formation rate ${J}_{M\varphi}$ as well [30], [33], [35]:(21)${J}_{M\varphi}={k}_{a}^{M}{\varphi}^{N}-{k}_{d}^{M}M\text{.}$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 $\frac{{k}_{a}^{M}}{{k}_{d}^{M}}=\frac{M}{{\varphi}^{N}}$. Now, if $\mathrm{\Phi}(=\mathit{NM}+\varphi )$ is the total concentration of surfactant monomers in the bulk, substitution of the initial total concentration ${\mathrm{\Phi}}_{0}$ results in:(22)$\frac{{k}_{a}^{M}}{{k}_{d}^{M}}=\frac{{\mathrm{\Phi}}_{0}-\mathit{CMC}}{N{(\mathit{CMC})}^{N}}\text{.}$Here, it should be noted that in equilibrium $\varphi =\mathit{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.

### 2.5. Evaporation

The evaporative volume flux ${w}_{e}$ is given by:(23)${w}_{e}=\frac{{D}_{v}{M}_{l}{p}_{\mathit{sat},l}}{\rho {R}_{u}T}\frac{\partial {\widehat{p}}_{l}}{\partial n}\sqrt{1+{\left(\frac{\partial h}{\partial r}\right)}^{2}}\text{.}$Here, ${D}_{v}$ is the vapor diffusivity in air, ${M}_{l}$ the molar mass of the liquid, ${p}_{\mathit{sat},l}$ the liquid saturation pressure, which is assumed constant, and ${\widehat{p}}_{l}=\frac{{p}_{l}}{{p}_{\mathit{sat},l}}$, with ${p}_{l}$ the local vapor pressure. The derivative $\frac{\partial}{\partial n}$ is the spatial derivative along the normal vector (pointing outwards) [21].

In order to calculate $\frac{\partial {\widehat{p}}_{l}}{\partial 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:(24)${\nabla}^{2}{\widehat{p}}_{l}=0\text{.}$The corresponding boundary conditions are given by:(25)${\widehat{p}}_{l}{|}_{z=h}=1\phantom{\rule{2em}{0ex}}\mathrm{for}\phantom{\rule{0.5em}{0ex}}r<R\text{;}$(26)$\frac{\partial {\widehat{p}}_{l}}{\partial z}{|}_{z=0}=0\phantom{\rule{1em}{0ex}}\mathrm{for}\phantom{\rule{0.5em}{0ex}}r>R\text{;}$(27)$\phantom{\rule{2em}{0ex}}{\widehat{p}}_{l}=\mathit{RH}\phantom{\rule{2em}{0ex}}\mathrm{for}\phantom{\rule{0.5em}{0ex}}(r,z)\to \infty \text{.}$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.

### 2.6. 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.

Name | Symbol | Definition | Expression |
---|---|---|---|

Evaporation number | $\mathrm{Ev}$ | $\frac{\mathrm{Evaporation}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}$ | $\frac{2}{\pi}\frac{\mu {D}_{v}{p}_{\mathit{sat},l}{M}_{l}}{H{({R}_{u}T)}^{2}{\mathrm{\Gamma}}_{0}\rho}(1-\mathit{RH})$ |

Transport number | $\mathrm{Tr}$ | $\frac{\mathrm{Desorption}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}$ | $\frac{{k}_{d}^{\mathrm{\Gamma}}{R}^{2}\mu}{{R}_{u}T{\mathrm{\Gamma}}_{0}H}$ |

Desorption number | $\mathrm{De}$ | $\frac{\mathrm{Desorption}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}{\mathrm{Adsorption}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}$ | $\frac{{k}_{d}^{\mathrm{\Gamma}}R}{{k}_{a}^{\mathrm{\Gamma}}}$ |

Surface Péclet number | ${\mathrm{Pe}}_{s}$ | $\frac{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}{\mathrm{Surface}\phantom{\rule{0.35em}{0ex}}\mathrm{diffusion}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}$ | $\frac{{\mathrm{\Gamma}}_{0}{R}_{u}\mathit{TH}}{{D}_{\mathrm{\Gamma}}\mu}$ |

Bulk Péclet number | ${\mathrm{Pe}}_{b}$ | $\frac{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}{\mathrm{Bulk}\phantom{\rule{0.35em}{0ex}}\mathrm{diffusion}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}$ | $\frac{{\mathrm{\Gamma}}_{0}{R}_{u}\mathit{TH}}{{D}_{\varphi}\mu}$ |

Micelle transport number | ${\mathrm{Tr}}_{M}$ | $\frac{\mathrm{Micelle}\phantom{\rule{0.35em}{0ex}}\mathrm{decomposition}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}$ | $\frac{{k}_{d}^{M}{R}^{2}\mu}{{R}_{u}T{\mathrm{\Gamma}}_{0}H}$ |

Micelle decomposition number | ${\mathrm{De}}_{M}$ | $\frac{\mathrm{Micelle}\phantom{\rule{0.35em}{0ex}}\mathrm{decomposition}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}{\mathrm{Micelle}\phantom{\rule{0.35em}{0ex}}\mathrm{formation}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}$ | $\frac{{k}_{d}^{M}{(\mathit{CMC})}^{1-N}}{{k}_{a}^{M}}$ |

Micelle Péclet number | ${\mathrm{Pe}}_{M}$ | $\frac{\mathrm{Marangoni}\phantom{\rule{0.35em}{0ex}}\mathrm{velocity}}{\mathrm{Micelle}\phantom{\rule{0.35em}{0ex}}\mathrm{diffusion}\phantom{\rule{0.35em}{0ex}}\mathrm{rate}}$ | $\frac{{\mathrm{\Gamma}}_{0}{R}_{u}\mathit{TH}}{{D}_{M}\mu}$ |

One may note that the traditional capillarity number $\mathit{Ca}=\mu {w}_{e}/{\sigma}_{0}$ is not listed in Table 1. This dimensionless number is omitted, because for any physically relevant case $\mathit{Ca}\ll 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.

## 3. 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.

### 3.1. 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 ${\theta}_{0}$ is always set to ${25}^{\circ}$ 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 coffee-ring 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 $\omega $ over the middle area where the Marangoni vortex tends to appear (around $R/2<r<5/6R$) shows that $\omega \u2a7d-1.2\xb7{10}^{-4}$ corresponds to the Marangoni regime, $-1.2\xb7{10}^{-4}<\omega \u2a7d-0.8\xb7{10}^{-4}$ to the transition regime and $\omega >-0.8\xb7{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 $\mathrm{Tr}<1$, both $1/\mathrm{Ev}$ and $\mathrm{Tr}$ have very similar effects on the occurence of Marangoni flow. This makes sense, because as $1/\mathrm{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 $\mathrm{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 $\mathrm{Tr}$ tends to promote Marangoni circulation is also found by Jung et al., who used a lattice-gas model [39].

For $\mathrm{Tr}>1$, the occurrence of Marangoni circulation becomes mostly limited by $1/\mathrm{Ev}$. Only if a certain threshold value of $1/\mathrm{Ev}$ is exceeded, Marangoni vortices can be formed. At the same time, $\mathrm{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/\mathrm{Ev}$ still is a limiting factor however, is not surprising. Decreasing surfactant concentration ${\mathrm{\Gamma}}_{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 $\mathrm{Tr}$ involves not only the desorption rate, but also the adsorption rate. If $\mathrm{Tr}$ is varied by changing ${k}_{d}$ the value of ${k}_{a}$ changes accordingly, given that $\mathrm{De}$ should remain constant. Increasing $\mathrm{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, $\mathrm{Tr}$ can be modified by changing ${\mathrm{\Gamma}}_{0}$ and by changing ${k}_{d}$. An increase in ${\mathrm{\Gamma}}_{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 $\mathrm{De}$ is shown. An increase of $\mathrm{Tr}$ still tends to promote Marangoni circulation, but this only holds approximately for $\mathrm{De}>10$. As the value of $\mathrm{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 $\mathrm{De}$ effectively is the solubility of the surfactant. As $\mathrm{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 $\mathrm{Tr}$ and low solubility to a low $\mathrm{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 $\mathrm{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 $\mathrm{Tr}$ and the high solubility to a high $\mathrm{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.

### 3.2. 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 ${\mathrm{Pe}}_{s}$ decreases to and becomes smaller than unity.

Similar behavior can be distinguished when the bulk Péclet number is varied. As ${\mathrm{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 ${\mathrm{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.

### 3.3. 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 ${\mathrm{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 ${\mathrm{De}}_{M}<1$ the value of ${\mathrm{De}}_{M}$ is largely irrelevant. The transition from coffee-ring regime to Marangoni regime is at that point only dependent on $\mathrm{Tr}$. For ${\mathrm{De}}_{M}>1$ however, a shift of the critical $\mathrm{Tr}$ to the left is observed. This shift can be explained by reformulating the definition of ${\mathrm{De}}_{M}$. Substitution of Eq. (22) in the definition of ${\mathrm{De}}_{M}$ yields:(28)${\mathit{De}}_{M}=\frac{{k}_{d}^{M}{(\mathit{CMC})}^{1-N}}{{k}_{a}^{M}}=\frac{N\xb7\mathit{CMC}}{{c}_{0}+{M}_{0}N-\mathit{CMC}}=\frac{\mathit{CMC}}{{M}_{0}}\text{.}$Here, it is used that ${c}_{0}=\mathit{CMC}$. Given this reformulation it becomes clear that varying the initial micelle concentration ${M}_{0}$ has a proportional effect on ${\mathrm{De}}_{M}$ and vice versa. A high value of ${\mathrm{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 ${\mathrm{De}}_{M}>1$.

In Fig. 4c it can be seen that the micelle Péclet number ${\mathrm{Pe}}_{M}$ has a similar effect on the flow dynamics as ${\mathrm{Pe}}_{s}$ and ${\mathrm{Pe}}_{b}$. As ${\mathrm{Pe}}_{M}$ is reduced, the transition from coffee-ring regime to Marangoni regime shifts to a higher $\mathrm{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].

## 4. 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–diffusion-adsorption 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 ($\mathrm{Tr},\mathrm{Ev},\mathrm{De}$) were analyzed for cases without diffusion and two additional parameters (${\mathrm{Pe}}_{s},{\mathrm{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 (${\mathrm{Tr}}_{M},{\mathrm{De}}_{M},{\mathrm{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.

Most preceding numerical studies focused either on the evaporation of droplets without surfactants [12], [21], [26], [27], [28] or on the evolution of nonvolatile droplets with surfactants [6], [29], [30], [31], [32], [33], [34], [35], [36], [37]. Only a few previous numerical studies involved both evaporation and surfactants [23], [38], [39]. However, this study is the first to numerically analyze the evolution of evaporating droplets with soluble surfactants and the corresponding internal flow patterns. Furthermore, various experimental studies consider the evaporation of droplets with surfactants [8], [17], [40], [41], [42], [43], [44], but actual flow visualisations are rare [24], [45]. Numerical work, like this study, can therefore be an attractive alternative.

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.

## CRediT authorship contribution statement

**R.T. van Gaalen:** Conceptualization, Methodology, Software, Investigation, Writing - original draft. **C. Diddens:** Methodology, Software, Supervision. **H.M.A. Wijshoff:** Funding acquisition, Supervision. **J.G.M. Kuerten:** Conceptualization, Writing - review & editing, Supervision.

## 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.

## Acknowledgements

This work is part of an Industrial Partnership Programme of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). This research programme is cofinanced by Canon Production Printing, University of Twente, Eindhoven University of Technology, and the “Topconsortia voor Kennis en Innovatie (TKI)” allowance from the Ministry of Economic Affairs.

## Supplementary material

The following are the Supplementary data to this article:

## References

- [1]Spray cooling heat transfer: the state of the artInt. J. Heat Fluid Flow, 29 (4) (2007), pp. 753-767, 10.1016/j.ijheatfluidflow.2006.09.003
- [2]Marangoni drying: a new extremely clean drying processLangmuir, 6 (11) (1990), pp. 1701-1703, 10.1021/la00101a014
- [3]High-speed and drop-on-demand printing with a pulsed electrohydrodynamic jetJ. Micromech. Microeng., 20 (2010), p. 095026, 10.1088/0960-1317/20/9/095026
- [4]Control of colloidal particle deposit patterns within picoliter droplets ejected by ink-jet printingLangmuir, 22 (8) (2006), pp. 3506-3513, 10.1021/la053450j
- [5]Inkjet printing-process and its applicationsAdv. Mater., 22 (6) (2010), pp. 673-685, 10.1002/adma.200901141
- [6]Recent advances in computational fluid dynamics relevant to the modelling of pesticide flow on leaf surfacesPest Manag. Sci., 66 (2010), pp. 2-9, 10.1002/ps.1824
- [7]Physical principles of pesticide behaviour, Vols. 1 and 2, Academic Press, London/New York (1980)
- [8]Droplet evaporation study applied to DNA chip manufacturingLangmuir, 21 (20) (2005), pp. 9130-9136, 10.1021/la050764y
- [9]Capillary flow as the cause of ring stains from dried liquid dropsNature, 389 (6653) (1997), pp. 827-829, 10.1038/39827
- [10]Pattern formation in drying dropsPhys. Rev. E, 61 (1) (2000), pp. 475-485, 10.1103/PhysRevE.61.475
- [11]Contact line deposits in an evaporating dropPhys. Rev. E, 62 (1) (2000), pp. 756-765, 10.1103/physreve.62.756
- [12]Numerical simulation of the drying of inkjet-printed dropletsJ. Colloid Interface Sci., 392 (2013), pp. 388-395, 10.1016/j.jcis.2012.09.063
- [13]Evaporating droplets on oil-wetted surfaces: suppression of the coffee-stain effectPNAS, 117 (29) (2020), pp. 16756-16763, 10.1073/pnas.2006153117
- [14]Suppression of the coffee-ring effect by shape-dependent capillary interactionsNature, 476 (2011), pp. 308-311, 10.1038/nature10344
- [15]Marangoni effect reverses coffee-ring depositionsJ. Phys. Chem. B, 110 (14) (2006), pp. 7090-7094, 10.1021/jp0609232
- [16]Overcoming the ‘coffee-stain’ effect by compositional Marangoni-flow-assisted drop-dryingJ. Phys. Chem. B, 116 (22) (2012), pp. 6536-6542, 10.1021/jp3009628
- [17]Surfactant-induced Marangoni eddies alter the coffee-rings of evaporating colloidal dropsLangmuir, 28 (11) (2012), pp. 4984-4988, 10.1021/la204928m
- [18]Buoyancy and Marangoni effects in an evaporating dropJ. Thermophys. Heat Transfer, 16 (2002), pp. 562-574, 10.2514/2.6716
- [19]On the effect of Marangoni flow on evaporation rates of heated water dropsLangmuir, 24 (17) (2008), pp. 9207-9210, 10.1021/la801294x
- [20]Marangoni convection instability in a sessile droplet with low volatility on heated substrateInt. J. Heat Mass Transfer, 117 (2017), pp. 274-286, 10.1016/j.ijthermalsci.2017.04.007
- [21]Modeling the evaporation of sessile multi-component dropletsJ. Colloid Interface Sci., 487 (2017), pp. 426-436, 10.1016/j.jcis.2016.10.030
- [22]Evaporation-triggered microdroplet nucleation and the four life phases of an evaporating Ouzo dropProc. Natl. Acad. Sci. USA, 113 (31) (2016), p. 8642, 10.1073/pnas.1602260113
- [23]Evaporation of sessile droplets laden with particles and insoluble surfactantsLangmuir, 32 (2016), pp. 6871-6881, 10.1021/acs.langmuir.6b01042
- [24]Surfactant-driven flow transitions in evaporating dropletsSoft Matter, 12 (2016), pp. 1593-1600, 10.1039/C5SM02354H
- [25]Control of stain geometry by drop evaporation of surfactant containing dispersionsAdv. Colloid Interface Sci., 222 (2015), pp. 275-290, 10.1016/j.cis.2014.08.004
- [26]Evaporating pure, binary and ternary droplets: thermal effects and axial symmetry, breakingJ. Fluid Mech., 823 (2017), pp. 470-497, 10.1017/jfm.2017.312
- [27]Detailed finite element method modeling of evaporating multi-component dropletsJ. Comput. Phys., 340 (2017), pp. 670-687, 10.1016/j.jcp.2017.03.049
- [28]Marangoni contraction of evaporating sessile droplets of binary mixturesLangmuir, 33 (2017), pp. 4682-4687, 10.1021/acs.langmuir.7b00740
- [29]Numerical modeling of drop coalescence in the presence of soluble surfactantsJ. Comput. Appl. Math., 293 (2016), pp. 7-19, 10.1016/j.cam.2015.04.013
- [30]Surfactant-enhanced rapid spreading of drops on solid surfacesLangmuir, 25 (24) (2009), pp. 14174-14181, 10.1021/la9019469
- [31]Surfactant-assisted spreading of a liquid drop on a smooth solid surfaceJ. Colloid Interface Sci., 287 (2005), pp. 233-248, 10.1016/j.jcis.2005.01.086
- [32]Effects of surfactant on droplet spreadingPhys. Fluids, 16 (3070) (2004), 10.1063/1.1764827
- [33]Surfactant-induced fingering phenomena beyond the critical micelle concentrationJ. Fluid Mech., 564 (2006), pp. 105-138, 10.1017/S0022112006001352
- [34]The spreading and stability of a surfactant-laden drop on a prewetted substrateJ. Fluid Mech., 554 (2006), pp. 5-24, 10.1017/S0022112005008104
- [35]On surfactant-enhanced spreading and superspreading of liquid drops on solid surfacesJ. Fluid Mech., 670 (2011), pp. 5-37, 10.1017/S0022112010005495
- [36]The effect of surfactant convection and diffusion on the evolution of an axisymmetric pendant dropletPhys. Fluids, 24 (2012), p. 062104, 10.1063/1.4729449
- [37]Film drainage between two surfactant-coated drops colliding at constant approach velocityJ. Colloid Interface Sci., 257 (2003), pp. 93-107, 10.1016/S0021-9797(02)00033-4
- [38]The evaporation of surfactant-laden droplets: a comparison between contact line modelsJ. Colloid Interface Sci., 579 (2020), pp. 888-897, 10.1016/j.jcis.2020.06.099
- [39]Surfactant effects on droplet dynamics and deposition patterns: a lattice gas modelSoft Matter, 13 (2017), p. 6529, 10.1039/C7SM01224A
- [40]Diffusion-controlled evaporation of sodium dodecyl sulfate solution drops placed on a hydrophobic substrateJ. Colloid Interface Sci., 362 (2011), pp. 524-531, 10.1016/j.jcis.2011.06.060
- [41]The coupling between evaporation and adsorbed surfactant accumulation and its effect on the wetting and spreading behaviour of volatile drops on a hot surfaceJ. Petrol. Sci. Eng., 51 (3–4) (2006), pp. 238-252, 10.1016/j.petrol.2005.12.008
- [42]Evaporation of droplets of surfactant solutionsLangmuir, 29 (2013), pp. 10028-10036, 10.1021/la401578v
- [43]Influence of surfactants on an evaporating drop: fluorescence images and particle deposition patternsLangmuir, 19 (2003), pp. 8271-8279, 10.1021/la030049t
- [44]Evaporation rate and development of wetted area of water droplets with and without surfactant at different locations on waxy leaf surfacesBiosyst. Eng., 106 (2010), pp. 58-67, 10.1016/j.biosystemseng.2010.02.004
- [45]Auto-production of biosurfactants reverses the coffee ring effect in a bacterial systemNat. Commun., 4 (2013), p. 1757, 10.1038/ncomms2746
- [46]Viscous-flow down a slope in the vicinity of a contact linePhys. Fluids A, 3 (1991), p. 515, 10.1063/1.858113
- [47]Free surface Stokes flow over topographyPhys. Fluids, 13 (2001), p. 2751, 10.1063/1.1401812
- [48]Two-dimensional simulations of flow near a cavity and a flexible solid boundaryPhys. Fluids, 18 (2006), p. 063103, 10.1063/1.2204061
- [49]On application of lubrication approximations to nonunidirectional coating flows with clean and surfactant interfacesPhys. Fluids, 22 (2010), p. 092102, 10.1063/1.3484276
- [50]Analysis of the microfluid flow in an evaporating sessile dropletLangmuir, 21 (9) (2005), pp. 3963-3971, 10.1021/la047528s
- [51]Analysis of the effects of Marangoni stresses on the microflow in an evaporating sessile dropletLangmuir, 21 (9) (2005), pp. 3972-3980, 10.1021/la0475270
- [52]Capillarity and Wetting PhenomenaSpringer, New York (2004)
- [53]Self-similar flow and contact line geometry at the rear of cornered dropsPhys. Fluids, 17 (072101) (2005), 10.1063/1.1946607
- [54]Density-driven flows in evaporating binary liquid dropletsPhys. Rev. Lett., 121 (2018), p. 184501, 10.1103/PhysRevLett.121.184501
- [55]Gravitational effect in evaporating binary microdropletsPhys. Rev. Lett., 122 (2019), p. 114501, 10.1103/PhysRevLett.122.114501
- [56]Thermodynamically consistent description of the hydrodynamics of free surfaces covered by insoluble surfactants of high concentrationPhys. Fluids, 24 (102107) (2012), 10.1063/1.4758476
- [57]Paradoxical coffee-stain effect driven by the Marangoni flow observed on oil-infused surfacesColloids Surf. A, 522 (2017), pp. 355-360, 10.1016/j.colsurfa.2017.03.009
- [58]Simple theory of “stick-slip wetting hysteresisLangmuir, 11 (3) (1995), pp. 1041-1043, 10.1021/la00003a057
- [59]On the surfactant mass balance at a deforming fluid interfacePhys. Fluids, 8 (3203) (1996), 10.1063/1.869098
- [60]Adsorption dynamics of surfactants at the air/water interface: a critical review of mathematical models, data, and mechanismsColloids Surfaces A Physicochem. Eng. Asp., 100 (1995), pp. 1-45, 10.1016/0927-7757(94)03061-4
- [61]the spreading of heat or soluble surfactant along a thin filmPhys. Fluids A, 5 (1993), pp. 58-68, 10.1063/1.858789
- [62]Gradient dynamics models for liquid films with soluble surfactantPhys. Rev. Fluids, 1 (2016), p. 083903, 10.1103/PhysRevFluids.1.083903
- [63]Foundations of Colloid ScienceOxford University Press, Oxford (1991)
- [64]Kinetics of adsorption from micellar solutionsAdv. Colloid Interface Sci., 95 (2002), pp. 237-293, 10.1016/S0001-8686(00)00085-3
- [65]Evaporation of a sessile droplet on a substrateJ. Phys. Chem. B, 106 (6) (2002), pp. 1334-1344, 10.1021/jp0118322
- [66]Surface instabilities and patterning at liquid/vapor interfaces: Exemplifications of the “hairy ball theoremColloids Interface Sci. Commun., 5 (2015), pp. 5-7, 10.1016/j.colcom.2015.04.003
- [67]Effects of Tween 20 and Tween 80 on the stability of Albutropin during agitationJ. Pharm. Sci., 94 (6) (2005), pp. 1368-1381, 10.1002/jps.20365
- [68]Critical Micelle Concentration of Aqueous Surfactant Systems, NSRDS-NBS 36US. Government Printing Office, Washington, DC (1971)
- [69]Evaporation of dilute sodium dodecyl sulfate droplets on a hydrophobic substrateLangmuir, 35 (32) (2019), pp. 10453-10460, 10.1021/acs.langmuir.9b00824
- [70]The yield stress – a review or ‘
*παντα**ρ∊ι*’ – everything flows?J. Non-Newtonian Fluid Mech., 81 (1999), pp. 133-178, 10.1016/S0377-0257(98)00094-9 - [71]A mechanism for non-Newtonian flow in suspensions of rigid spheresJ. Rheol., 3 (1959), pp. 137-152, 10.1038/10.1122/1.548848
- [72]Intermolecular and Surface ForcesElsevier inc., London (2011)
- [73]Fabricating colloidal particles with photolithography and their interactionsat an air-water interfacePhys. Rev. E, 62 (2000), p. 951, 10.1103/PhysRevE.62.951