The evaporation of surfactant-laden droplets: A comparison between contact line models Journal of Colloid and Interface Science

Hypothesis: There are two different sharp-interface models for moving contact lines: slip models and pre- cursor ﬁlm models. While both models predict a mostly constant contact angle during the evaporation of pure droplets, it is expected that they behave differently when surfactants are present, because of the inherent dissimilarities in their respective interface deﬁnitions. Simulations: Both contact line models are numerically implemented using lubrication theory to analyze evaporating droplets. A convection-diffusion equation is implemented for insoluble surfactants. For pure droplets the models are compared with experiments performed by Nguyen et al. (2012). Findings: The two contact line models show results comparable to the experiments with pure droplets. If insoluble surfactants are present, the slip model increasingly shows pinning-like behavior as the initial surfactant concentration is increased. This ‘quasi-pinning’ is found to be consistent with experimental results in literature. The precursor ﬁlm model, in contrast, shows no signiﬁcant change when surfactants are added. This lack of change is a result of surfactant ﬂowing from the droplet into the precursor ﬁlm and vice versa. Whilesuggesting potentialsolutions tothis unphysicalbehavior,


Introduction
Sessile droplets are a common phenomenon. It can rightly be said that these droplets are encountered in a broad spectrum of contexts, ranging from everyday situations, as raindrops on a window and coffee spilled over a table, to high-tech applications, as inkjet printing [1] and spray cooling [2]. Especially the technological significance of sessile droplets has aroused a widespread interest in the subject. Increasingly detailed experiments and models have allowed for a much better understanding of the physical behavior of droplets in many configurations and on different scales. However, several challenges and questions still remain, which require continued research.
One topic that is still actively investigated is the drying of sessile droplets. This drying process is typically controlled by the diffusion of vapor from the droplet interface to the environment, as was demonstrated by Deegan et al. [3,4] and Popov [5]. As a result, a nonuniform evaporation rate is induced, with a singularity at the contact line if the contact angle is small. The wetting behavior of sessile droplets can be divided into partial wetting, for which the droplet has an associated equilibrium contact angle [6], and complete wetting, for which no equilibrium contact angle exists, which results in indefinite spreading of the droplet [7]. Grounded on these principles, numerous analytical and numerical studies have been carried out [8][9][10].
During evaporation, the contact line typically shows two different types of behavior: either the contact angle decreases while the contact line remains pinned or the contact line moves while the contact angle remains constant. Also, a combination of these two types can occur in the form of stick-slip behavior [11]. Pinning of the contact line can occur if the contact angle is small or if the substrate is relatively rough. This is the case in many practical applications [12][13][14]. Constant contact angle evaporation on the other hand, typically occurs if the substrate is smooth or the contact angle is large. Both modes of evaporation can be modelled with hydrodynamic equations, but if the contact line moves an additional description of the nanometer scale around the contact line is required to resolve the contact line singularity [15].
One possible way of dealing with the contact line singularity is by assuming a slip length, that allows for some flow parallel to the substrate [16][17][18]. This is both experimentally and numerically a well established method, but a drawback is that it only introduces a length scale and not an energy scale that describes the interaction with the substrate. Practically, this means that an additional boundary condition needs to be introduced for the microscopic contact angle to close the problem. This boundary condition is usually an equation that couples the contact line velocity to the contact angle. For more information on this subject, the reader is referred to the review by Snoeijer and Andreotti [19].
Another method of resolving the singularity is by introducing a precursor film that covers the substrate around the droplet. In that case, the contact line is defined on the free surface of the film, where slip is allowed, thus resolving the contact line singularity [6,20,21]. In order to keep the precursor film stable, a disjoining pressure is introduced, which accounts for the molecular interactions of the fluid with the substrate [22][23][24][25]. An advantage of the precursor film model is that it effectively resolves the contact line singularity without the need for additional boundary conditions. A disadvantage however, is that precursor films are rarely encountered under partial wetting conditions [19,26,27].
Given these two contact line models, it is considered of interest to compare them in the context of evaporating droplets to see whether both models give rise to similar results. While some studies exist that compare these two models (e.g. [28][29][30] in the context of spreading and [31] for a flow driven by a surface tension gradient), a comparison has never been made for evaporating dro-plets. Depending on the outcome, such a comparison is relevant since it guides researchers in making a substantiated decision with regard to the applicability of these models.
Furthermore, besides comparing the two contact line models for pure droplets, it is also relevant to compare the models with the inclusion of insoluble surfactants. Surfactants decrease the local surface tension of an interface. In principle, all surfactants are to some extent soluble, which implies that there exists a surfactant concentration both at the free surface and in the bulk of the liquid. Between these concentrations surfactant is transferred through adsorption and desorption. However, for many practically relevant surfactants the bulk solubility is very small (e.g. oleic acid in water [32,33]). For these surfactants the bulk concentration can be disregarded in which case one speaks of 'insoluble surfactants' [34].
While several numerical studies exist that consider pure, volatile droplet dynamics (e.g. [21,35]) or the dynamics of nonvolatile droplets and surfactants (e.g. [36][37][38]), numerical studies that analyze both the evolution of evaporating droplets and surfactants are less common. This is only done in the work of Karapetsas et al. [39], who model the evaporation of a complete wetting droplet with insoluble surfactants. Thus, it is of great relevance to analyze the evaporation of surfactant-laden droplets under partial wetting conditions, especially given the wide application of surfactants.
In the present work, a comparison is made between the precursor film model and slip length model in the context of evaporating sessile droplets with insoluble surfactants. For both pure droplets and droplets with insoluble surfactants the volume, radial and contact angle evolution calculated by both models is analyzed and compared. Also, both models are quantitatively compared with experiments for pure droplets (i.e. droplets without surfactants) and checked for consistency with experiments for droplets with insoluble surfactants.
It is hypothesized that for pure droplets both the precursor film model and slip model will yield similar results, because both involve evaporation with a mostly constant contact angle. Furthermore, it is expected that during the evaporation of surfactant-laden droplets the surfactant concentration will increase as a result of the decreasing interfacial area. This will decrease the contact angle as the surface tension decreases with increasing concentration. Because of the inherent dissimilarities in interface definitions between the models, it is foreseen that the degree with which the surfactant concentration increases differs per model. This article is organized as follows: first, in Section 2 the applied numerical method is presented, which is based on lubrication theory. The droplet height evolution equation and interfacial surfactant transport equation are introduced and an evaporation model is derived. Extra attention is paid to the formulation of the involved pressure terms and to the boundary conditions in the slip length model. After that, in Section 3 the results of the numerical simulations are shown and the effects of both models on the volume, radial and contact angle evolution are analyzed and compared. In the final section, conclusions are drawn based on the found results.

Mathematical model
In this section the mathematical equations are introduced that model the evaporating droplet. First, the evolution of the height profile is presented and described and after that the convectiondiffusion equation for the surfactants at the interface. Subsequently, a separate subsection is dedicated to the formulation of the pressure terms in the system, including the disjoining pressure that is related to the precursor film. Lastly, the slip model is outlined and the evaporation model is explained.

Lubrication equation
A sessile droplet is considered with mass density q, dynamic viscosity l and (liquid-gas) surface tension r. For this study, only droplets much smaller than the capillary length l r ¼ ðr=qgÞ 1=2 are regarded. Therefore, the influence of gravity g can safely be neglected [20,40]. Another consequence of the small length scales involved is a low Reynolds number. This results in a viscous flow where the inertia terms can be neglected. Lastly, cylindrical coordinates (r; /; z) are used with the assumption of axisymmetry ( @ @/ ¼ 0) and no swirl (U / ¼ 0) and, most importantly, 'flat' droplets are assumed, which allows one to express the evolution equation of drop height hðr; tÞ as a function of r: Here, p is the excess pressure in the droplet and w e is defined as the height change rate due to evaporation. Both these parameters will be quantified later on, in SubSection 2.3 and SubSection 2.5, respectively. In the evolution equation, a pressure driven term can be recognized and a term driven by the surface tension gradient, which accounts for the Marangoni effect [41]. For a derivation of Eq. (1), see Subsection S.1.1.
The boundary and initial conditions that hðr; tÞ is subjected to are given by: Here, R is the radius of the droplet and h 0 is the initial droplet profile, which is obtained by setting w e ¼ 0 and relaxing the droplet to an equilibrium solution.

The evolution equation of the local surfactant concentration
Cðr; tÞ can be expressed in radial coordinates as: Here, U t is the fluid velocity tangential to the liquid-air interface and D s the surface diffusion coefficient. Furthermore, @ @s denotes the partial derivative tangential to the interface in radial direction. The first term on the right-hand side accounts for convective transport tangential to the interface, the second term is a source-like term that takes into account deformation of the interface, the third term is a diffusion term and the last term compensates for the displacement of the surface coordinates that move along the normal of the surface. A derivation of Eq. (5) is given in Section S.2.
Surfactant concentration Cðr; tÞ is subject to the following boundary and initial conditions: Here, D denotes the domain boundary for the surfactant, which is where @ r h % 0 for the precursor film model and D ¼ R for the slip length model.
The surface tension equation of state used in this present work, is the dilute limit relation given by In this relation r C is the 'efficiency' of the surfactant which is defined as r C ¼ R u T with R u the universal gas constant and T the temperature.

Disjoining pressure and precursor film
The pressure p in the droplet is given by Derjaguin's equation [21,22,42]: where p L is the Laplace pressure and P the disjoining pressure. Both these pressure terms are independent of z, which is in accordance with Equation (S.4). The Laplace pressure is given by Contrary to the traditional formulation of the Laplace pressure, surface tension r is also part of the derivative in Eq. (11), rather than that it is written outside the brackets. This is because r depends on the local surfactant concentration and thus it depends on r. As shown by Thiele et al. [43] the surface tension is required to be part of the radial derivative in that case. Disjoining pressure P is in this model given by Here, B; n and m are positive constants (with n > m) and h s is the precursor film height for which P equals zero. The height h s does not necessarily correspond to the real precursor film height, but only represents the zero-pressure height. This formulation of the disjoining pressure consists of an attractive and a repulsive component (Bðh s =hÞ m and Bðh s =hÞ n respectively), corresponding to the attractive and repulsive molecular forces that govern the physics of the precursor film [25]. In literature, one often encounters disjoining pressure formulations that only include the attractive forces (see e.g. Karapetsas et al. [39]). In these cases the droplet will spread indefinitely, because there is no equilibrium contact angle corresponding to this formulation as shown by Starov and Velarde [6]. In this work, n and m are chosen to be 5 and 4 respectively. Higher values of n and m would require a finer grid without any significant changes in the drop profile, while lower values of n and m would cause significant deviations in the drop profile from the spherical cap shape. The value of B can be related to the equilibrium contact angle h e , which is elaborated in Section S.3.

Slip length model
Another way of solving the contact line singularity is by assuming a slip length (see e.g. [16]). The height evolution equation given by Eq. (1) then becomes: singularity, but still requires a complementary constitutive relation that prescribes the contact line velocity [19]. In this study, the contact line velocity is obtained by means of the relation: Here, h is the current contact angle, a a power-law index with a range of 1-3 and k a typical sensitivity of the contact line position to deviations of the contact angle from equilibrium with the dimension of velocity. Increasing a typically results in a less constant contact angle during evaporation, since h À h e is generally smaller than unity (the units are radians). An increase of k results in a more constant contact angle, since the radius will respond faster to any deviations from h e . In this work, a is set to 1 and k is set to 1:3 Â 10 À5 m=s. These values are empirically determined by comparison with experimental data from Nguyen et al. [44]. This data is also used for validation of the numerical models in SubSection 3.1. Note that for small deviations from the equilibrium contact angle, Eq. (14) with a ¼ 1 is in agreement with the Cox-Voinov law.
The current contact angle h can be calculated from the slope @ r h at the contact line and, in case of surfactants, h e can be computed by means of Young's equation: ð15Þ Here, h 0 is the initial contact angle.

Evaporation model
The evaporative mass flux of water _ m w is given by: with D wa the water vapor diffusivity in air [45], M w the molar mass, p sat;w the homogeneous saturation pressure, which is assumed to be constant, andp w ¼ p w p sat;w with p w the local vapor pressure [35].
Regarding the vapor field, Deegan et al., [4] and Hu and Larson [46] showed that the vapor diffusion for single component droplets can be considered as instantaneous. Combined with the assumption of no air flow, this results in the following Laplace equation: The corresponding boundary conditions of Eq. (18) are: Eq. (19) represents the saturated vapor at the liquid-air interface, Eq. (20) the prohibition of vapor penetration in the substrate and Eq. (21) the relative humidity in the surroundings (where RH w is the relative humidty). For the precursor film model, the radial coordinate corresponding to R is defined to be the point where h ¼ 1:2h s , which results in the droplet having approximately the same effective radius as if it would have been a spherical cap.
The evaporation velocity w e in Eq. (1) can subsequently be calculated by: Similarly to Eq. (11), the argument of the square root makes the expression a good approximation for higher slopes than it would when approximated as unity. An elaboration of the numerical method is given in Section S.4.

Results and discussion
In this section the numerical results are shown and discussed. First, both the precursor film model and slip length model are compared to experimental results, all for pure water droplets. After that, the effects of insoluble surfactants on the results of these two cases are investigated and explained.

Pure droplets
First, the case of pure droplets is investigated. Simulations are carried out for 1.45 microliter water droplets. Furthermore, all physical properties of water are taken at 25 C and the relative humidity RH w ¼ 0:55. These parameters are chosen in order to compare the simulations with the experiments that were performed by Nguyen et al. [44]. They carried out experiments for water droplets on an atomically smooth Oct-Silicon substrate and published data both in the pinned contact line regime and constant contact angle regime. Precursor film height h s is set equal to the height of one grid cell (%9.4 lm) and slip length b ¼ 100 nm.
Although physical values of these two parameters tend to be smaller, simulations and literature show that any smaller values result in negligible changes in the solution, but require a finer grid [21].
The initial contact angle was set at h 0 ¼ 35:0 for both models. This yielded nearly identical height profiles as can be seen in Fig. 1a (h 0 ¼ 31:5 will be discussed later). For the precursor film model, the contact angle is defined as the contact angle of a spherical cap with the same height and volume (after subtraction of the volume of the precursor film). Calculating the contact angle directly from the height profile, for example from the highest slope, would result in an underestimation of the effective contact angle, since the drop blends into the precursor film around R.
As can be seen in Fig. 2 both the precursor film model and the slip length model agree rather well with each other and with the experiments. The differences between both models are minute in terms of volumetric evolution and radial evolution and the dissimilarities that can be observed are likely caused by the fact that both models give slightly different drop shapes during evaporation. In Fig. 1b it is shown that at t ¼ 600 s the slip length model yields a flatter, lower-contact-angle drop than the precursor film model. This behavior can also be observed in Fig. 2, especially in the contact angle evolution. The precursor film model predicts a nearly constant contact angle that only decreases at the end of the drying process, while the slip length model also predicts an initial decrease. This initial decrease can be explained by Eq. (14), which shows that the change in R is proportional to the deviation from h e . As a result, the contact angle will keep decreasing due to evaporation, until @ t R is large enough to maintain a quasi-steady, receding contact angle. This initial decrease in contact angle is also observed in the experiments.
Compared with the experiments, two discrepancies can be distinguished in the volume evolution. Firstly, at the start the experimental results show a slight deviation from the model results. This deviation can be attributed to the transition between the pinned contact line regime, which occurs in the experiments before the initial time of the simulations and is not shown in the figure, and the moving contact line regime, which is actually shown in the figure. Because the radius is nearly constant during the transition, the evaporation occurs faster than it would for a moving contact line. The slip length model partially captures this effect, as can be seen in Fig. 2c, but the precursor film model has a nearly constant contact angle until the end of the drying process. Therefore, more comparable results would be obtained for the precursor film model if the initial contact angle is set equal to the receding contact angle that is observed in the experiments (% 31:5 ). Therefore, it is decided to use h 0 ¼ 31:5 for the precursor film model to obtain more representative and comparable results. This will cause the drop profiles to deviate in the beginning of the evaporation process, as can be seen in Fig. 1a, but they will be more similar at later stages, as illustrated in Fig. 1b.
The second discrepancy between the models and experiments can be seen at the end of the drying process (around t ¼ 1000 s), where the experiments again show nonlinear behavior. This  nonlinear behavior can also be explained by the decrease of the contact angle, which accelerates the drying process. Both models show qualitatively the same behavior, but predict that the contact angle is constant for a longer duration. Unfortunately, how large the difference in contact angle between the models and the experiments is after t ¼ 1000 s cannot be said with certainty, since Nguyen et al. did not report any data for that time (probably because the contact angle measurements became increasingly less accurate when h ! 0).
In summary, it is concluded that both the precursor film model and the slip length model match the experiments reasonably, slightly favoring the slip model when it comes to contact angle evolution. Of course, it is duly noted that the agreement is partly a result of fitting moving contact line parameters (n and m for the precursor film model, a and k for the slip model) to the experiments, as reported in SubSections 2.3 and 2.4.

Droplets with insoluble surfactants
In the previous subsection a pure water droplet was considered for a precursor film model and a slip length model. However, many applications involve surfactants, either intentionally or in the form of impurities. Therefore, in this subsection sessile droplets with insoluble surfactants are studied and again the two contact line models are compared.
Surfactant efficiency r C ¼ R u T is set to 2:48 kJ/mol, which corresponds to a temperature of 25 C, and surface diffusion coefficient D s is set to 4:33 Â 10 À10 m 2 =s, which corresponds to 2D Einsteinian diffusion of a molecule with an effective radius of five times the Van der Waals radius of a carbon atom (850 pm) [47]. All simulated droplets have the same initial equilibrium contact angle (35:0 for the slip length model and 31:5 for the precursor film model), regardless of the initial surfactant concentration. In Fig. 3, it can be observed that for the slip length model the addition of surfactants has a significant influence on the drop evolution. Higher concentrations result in a shorter drying time (see Fig. 3a), which is caused by the fact that the equilibrium contact angle decreases as the surfactant concentration increases due to the shrinking interface. This also results in apparent contact line pinning, as can be seen in the radial evolution: any evaporated volume results in an increased surfactant concentration, which subsequently reduces the contact angle accordingly, while keeping the radius nearly constant. Thus, it can be seen that for higher initial surfactant concentrations the contact angle decreases at a higher rate, while the radius tends to be more constant. This also causes a relatively higher concentration increase for the lower initial concentrations, since there the interfacial area shrinks faster due to the decreasing radius.
This 'quasi-pinning' behavior is consistent with results found in literature. As reported by Truskett and Stebe [48], pinning behavior is strongly promoted with the addition of surfactants. While this increased pinning can also be attributed to the fact that surfactants decrease the initial contact angle, which promotes pinning as well, the quasi-pinning observed in the simulations can be an additional relevant mechanism. Similar results have been reported for soluble surfactants by Dugas et al. [49], who show improved pinning by surfactants, and by Sefiane [50], who, besides increased pinning, also shows a decrease in contact line retraction rate with the addition of surfactants. Both these results are also consistent with the quasi-pinning behavior observed in the simulations, because the mechanism is essentially the same for soluble surfactants: a reduction in available space for surfactants results in a concentration increase, which subsequently results in a further reduced equilibrium contact angle.
The occurrence of quasi-pinning is also predicted by Semenov et al. [51]. They show both experimentally and analytically that the volume evolution of evaporating, surfactant-laden droplets scales as V 2=3 $ t for contact angles higher than approximately 45°, and predict that it accelerates for contact angles lower than 45°.
As a side note, it is not taken into account here that surfactants might also inhibit evaporation (e.g. as observed by [52], while not observed by [48]). However, this would only affect the time scale. The observed quasi-pinning of the contact line will not be affected by including a reduced evaporation rate, since the pinning is a result of the equilibrium contact angle being reduced by the increased surfactant concentration.
Another relevant observation is that in none of the simulations the typical circulating flow occurs, that is often caused by the Marangoni effect (e.g. see Fig. 4a. Note that it is possible to model flow circulations with lubrication theory despite the small contact angle assumption [53].). In first instance, this may seem counterintuitive given the frequency with which surfactant-induced Marangoni vortices are reported in literature (e.g. see [54]). However, we will show that, in this case, the occurrence of a Marangoni vortex is physically impossible, regardless of the initial surfactant concentration.
The argument is the following. Given that evaporation is strongest at the contact line, there will be a higher curvature at the apex compared to the contact line. The corresponding pressure gradient causes the typical outward, 'coffee-ring' flow [55,56]. Furthermore, because of the same curvature effect, the surfactant concentration increases most rapidly at the drop apex. This source-like effect is a direct result of the second term on the right-hand side of Eq. (5). The movement of the contact line does not change this effect, as implied by Equation (S.25). Thus, the surfactant concentration at the apex can never be smaller than at the contact line due to shrinkage of the interface, although a positive concentration gradient towards the contact line is required for a circulating flow.
Thus, for Marangoni vortices to occur the required surface tension gradient should be caused by convective transport effects. If any surfactant is transported towards the contact line by the coffee-ring flow, a negative surface tension gradient could arise. However, this surface tension gradient can never be sufficient to cause a circulating flow, because as soon as it would be large enough for a negative interfacial velocity, it would cancel the outward flow at the interface. Subsequently, without an outward flow, the concentration at the contact line will not increase anymore. Thus, it can be concluded that the surfactant tends to reduce the interfacial flow, but it can never cause a Marangoni vortex. Fig. 4b is consistent with this conclusion: the negative pressure gradient induces an outward flow that is partly diminished due to the resulting concentration gradient.
For Marangoni vortices to occur, other physical effects need to be taken into account. For example, if the surfactant is soluble, it can also be transported through the bulk rather than only along the interface. Thus, circulations can then be caused by surfactant continuously flowing through the bulk towards the contact line, adsorbing at the interface, flowing along the interface towards the drop apex and desorbing into the bulk. Alternatively, if an additional driving force is involved, like spreading, circulations can occur with insoluble surfactants [39]. Lastly, temperature effects and other solutal effects can also cause Marangoni vortices [35].
When observing the drop evolution as predicted by the precursor film model (Fig. 5), a major flaw of this model is revealed. Addition of surfactants has effectively no effect on the volume and radial evolution and higher surfactant concentrations hardly reduce the contact angle during evaporation. This is caused by the fact that the droplet is not well defined in the precursor film model and therefore there is no mechanism to keep the surfactant from flowing into the precursor film. Therefore, while the droplet shrinks, surfactant is either left behind in the precursor film or flows back into the precursor film to reduce the concentration gradient. This backflow gives rise to the tiny concentration increase that can be seen in Fig. 5d. This is not an issue that can easily be solved within the context of the precursor film model. An obvious way would be to introduce an artificial no-flux boundary condition at the same point where the evaporation is cut off. However, this only solves the problem of surfactant flowing out of the droplet, but does not cause the surfactant to move with the contact line. A more physical method would be to involve a penalty function that assigns a relatively high energetic cost to surfactant in the precursor film, for example as was done by Thiele et al. [57]. A major disadvantage of this method, however, is the resulting severe numerical instability. This makes penalty functions more suitable for equilibrium cases. An alternative strategy could be to introduce a source function that creates surfactant as the effective drop radius decreases. This would in principle solve the problem of too low concentration inside the droplet, but is of course highly unphysical.
In summary, the slip length model yields reasonable results that can be understood intuitively and are consistent with literature, while the precursor film model has some serious defects when surfactants are involved. Possibly, in future work methods will be derived that use disjoining pressure without precursor films as suggested by Colinet and Rednikov [58]. This would indeed solve the backflow problem associated with precursor films, but until then the slip length model is preferred.

Conclusion
A numerical model was developed using lubrication theory to simulate the evolution of evaporating droplets [59][60][61][62]. Insoluble surfactants were included in the model by means of an interfacial convection-diffusion equation [43,[63][64][65]. In order to incorporate moving contact lines, both a precursor film model and a slip length model were implemented. For both models, the numerical results were compared to experimental results with pure liquids obtained by Nguyen et al. [44]. This showed reasonable agreement, slightly favoring the slip length model.
In case of insoluble surfactants, the slip length model showed that 'quasi-pinning' of the contact line can occur. This behavior is a direct result of the fact that a shrinking interface causes the surfactant concentration to increase, which subsequently reduces the equilibrium contact angle. Thus, at some point any evaporated volume results in a decreasing contact angle, effectively pinning the contact line. Quasi-pinning was found to be consistent with experimental results obtained by Truskett and Stebe [48], who show that insoluble surfactants promote pinning. Furthermore, similar experimental results were obtained for soluble surfactants by Dugas et al. [49] and by Sefiane [50]. The possibility of quasi-pinning was also inferred by Semenov et al. [51]. Besides quasi-pinning, it was also shown that surfactants tend to reduce the bulk flow close to the interface, but will never cause the typical circulating flow pattern that is often observed when surfactants are involved.
The precursor film model, on the contrary, showed no significant differences in drop behavior with the inclusion of insoluble surfactants. This is a result of surfactant being left behind in the precursor film as the contact line moves and even surfactant flowing back from the droplet into the film. Several methods to counteract this unphysical behavior were discussed, but none are considered desirable or physical.
These results agree with the hypothesis that was made, namely that both contact line models would produce similar results for pure droplets, but different results when surfactants are included. Furthermore, the contact angle indeed decreases over time as surfactants are included, especially for higher initial concentrations.
From the perspective of progress and innovation the relevance of the findings is twofold. First, from the comparison between contact line models it is concluded that in the context of pure, evaporating, sessile droplets both models perform nearly equal, slightly favoring the slip length model in comparison with experiments. Furthermore, with the inclusion of surfactants the slip length model is preferable over the precursor film model given the current state of the art. Also, the pinning-like behavior exhibited by the surfactant-laden drop is of potential use as a mechanism for controlling the final drop radius in technologies such as inkjet printing, pesticide sprays, the fabrication of DNA/RNA microarrays and many more. This may be possible by fine tuning the initial surfactant concentration with regard to the desired final radius.
Previous comparisons between contact line models only included nonvolatile droplet dynamics [28][29][30][31], so a numerical comparison for volatile droplets is a welcome addition to the existing literature, especially given the increasing interest in evaporating droplets with moving contact lines [21,35,66,67]. Preceding numerical studies about surfactant dynamics mostly involved nonvolatile drops [16,17,[36][37][38][68][69][70][71][72][73], while only one numerical study known to the authors considers both surfactant and volatile droplet dynamics (Karapetsas et al. [39]). Since the latter treats droplets under complete wetting conditions the results of the present work, which considers partial wetting conditions, can be seen as complementary.
Further research opportunities lie in the modeling of contact angle hysteresis. It is generally known that during evaporation, water droplets first tend to go through a pinned contact line regime until the contact angle equals a receding contact angle h r . Only then will the contact line start moving, while keeping a constant contact angle, until it passes through a third and final regime where both the contact radius and contact line decrease [44,74]. In this work, only the latter two regimes are considered, so it may be of interest to modify the contact line models to include the pinning regime as well. For the slip length model, this can be done by setting @R=@t ¼ 0 for h > h r and for the precursor film model, it might be possible to introduce a prefactor B that depends on h. An illustration of drop evaporation with contact line hysteresis is shown in Fig. 6. Surfactants could be included by making h r dependent on surfactant concentration C.
Regarding evaporating droplets with surfactants in general, a fitting next step is an extension of the model to soluble surfactants, micelles and adsorption on the substrate, similar to the work of Karapetsas et al. [71]. This will render the models increasingly rel-

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.