Pattern formation during the evaporation of a colloidal nanoliter drop: a numerical and experimental study

An efficient way to precisely pattern particles on solid surfaces is to dispense and evaporate colloidal drops, as for bioassays. The dried deposits often exhibit complex structures exemplified by the coffee ring pattern, where most particles have accumulated at the periphery of the deposit. In this work, the formation of deposits during the drying of nanoliter colloidal drops on a flat substrate is investigated numerically and experimentally. A finite-element numerical model is developed that solves the Navier–Stokes, heat and mass transport equations in a Lagrangian framework. The diffusion of vapor in the atmosphere is solved numerically, providing an exact boundary condition for the evaporative flux at the droplet–air interface. Laplace stresses and thermal Marangoni stresses are accounted for. The particle concentration is tracked by solving a continuum advection–diffusion equation. Wetting line motion and the interaction of the free surface of the drop with the growing deposit are modeled based on criteria on wetting angles. Numerical results for evaporation times and flow field are in very good agreement with published experimental and theoretical results. We also performed transient visualization experiments of water and isopropanol drops loaded with polystyrene microspheres evaporating on glass and polydimethylsiloxane substrates, respectively. Measured evaporation times, deposit shapes and sizes and flow fields are in very good agreement with the numerical results. Different flow patterns caused by the competition of Marangoni loops and radial flow are shown to determine the deposit shape to be either a ring-like pattern or a homogeneous bump.

In biology, the evaporation of colloidal drops is used for depositing or organizing biological materials such as proteins and DNA [1]- [7]. Colloidal deposition and crystallization [8]- [14] can also be used to manufacture micro-and nanowires [15,16], nanocrystals [17] and explosive crystalline layers [18]. Figure 1 shows that the patterns left by an evaporating colloidal drop can exhibit a ring-like structure [19] or more complex features such as a network of polygons [20], hexagonal arrays [17] or a uniform deposit [21,22]. The pattern variety echoes the complex, coupled and multi-scale transport phenomena occurring during the evaporation of a colloidal drop: fluid dynamics is transient and can be influenced by the impact and associated interfacial deformation by Marangoni stresses, by wetting and the shrinking free surface; heat transfer occurs by convection inside the drop and conduction in the substrate, with a latent heat contribution at the evaporating free surface; mass transfer occurs through diffusion of liquid vapor in the atmosphere, advection-diffusion of particles in the drop and short-range forces between the particles and the substrate. In 1997, Deegan et al explained the ring deposit formation [19,23,24] by a strong radial flow carrying particles toward the pinned wetting line, where evaporative flux is the highest due to the wedge geometry. Their theory [19,23,24] was based on the lubrication approximation. For thin drops with spherical cap shapes, they provided an analytical expression for the local evaporative flux. Maenosono et al [25] observed that the growth of a nanoparticle ring occurs in two stages, a ring buildup stage with a pinned contact line, followed by a stage where the wetting line recedes. They predicted the ring growth dynamics analytically and found a reasonable agreement with experiments in terms of ring growth rate and final width. Onoda and Somasundaran [26] showed that the chemical or physical heterogeneities of the substrate influence the final deposition pattern. Also, Marangoni convection inside the drop due to the variation of surface tension along the free surface can affect the deposit formation and the flow fields, as shown by Truskett and Stebe [20]. In 2006, Hu and Larson [22] observed that thermal Marangoni convection induces preferential deposition at the center of the droplet rather than at the periphery, for the case of microliter octane drops evaporating on glass coated with perfluorolauric acid (PFLA). Very recently, Ristenpart et al [27] showed analytically that the ratio between substrate and droplet thermal conductivities controls the direction of Marangoni convection inside an evaporating drop. They also showed experimentally that deposit patterns depend on the Marangoni flow direction. The formation of multiple rings has been investigated by Shmulyovich et al [28] during the drying of 1-70 µl water droplets containing 0.88 and 3.15 µm latex particles. The ring formation occurred as a succession of wetting line pinning, ring formation and depinning events. Deegan observed experimentally [23] that the time at which the wetting line of a colloidal drop depins (defined as depinning time) from the deposit corresponds to about 40-80% of the total evaporation time. In a very recent article, Zheng [29] showed theoretically that depinning occurs at about 50% of the total evaporation time, independently of the initial concentration of the particles in the drop.
Recently, numerical studies have been reported that study the evaporation of a sessile drop of pure liquid on a solid substrate. Hu and Larson [30,31] developed a finite-element model using the lubrication approximation for the evaporation of a sessile drop on a glass substrate with a pinned wetting line. They assumed that the free surface of the drop remains a spherical cap throughout the evaporation. In an extension of this model [32], the same authors reported Figure 1. A multiplicity of deposits can be obtained after the drying of a colloidal drop: (a) ring-like pattern from an aqueous drop containing 60 nm polystyrene spheres on titanium substrate [90] (with permission from ACS); (b) multiple rings from a µl water drop containing 1 µm polystyrene microspheres on glass (our work); (c) fingering at wetting line obtained from a µl isopropanol drop with 1 µm polystyrene microspheres on glass (our work); (d) uniform deposition pattern of 60 nm hydroxyapatite particles from an aqueous drop on a titanium disk [91] (with permission from ACS); (e) hexagonal cells from surfactant-laden aqueous drop containing polystyrene microspheres on hydrophobic octadecyltricholorosilane (OTS) substrate [20] (with permission from ACS). 6 that Marangoni convection changes the internal flow pattern from radially outward to counterclockwise. Also, a detailed numerical model for the evaporation of pure microliter water drops was presented by Ruiz and Black [33]. They solved the Navier-Stokes and energy equations with consideration of thermocapillary stresses. In this model, the wetting line was pinned and the evaporative flux was assumed to depend on the substrate temperature but not on the radial distance. Studies by Mollaret et al [34], Girard et al [35] and Widjaja et al [36] followed a similar approach assuming a pinned wetting line but solved the vapor diffusion equation to obtain an accurate evaporative flux.
Regarding colloidal drop evaporation, fewer numerical studies have been reported. Within the framework of the lubrication theory, Fischer [37] studied the advection of particles in an evaporating drop, and numerically found transient and local particle distributions within the droplet. Hu and Larson [22] computed the deposition of PMMA particles based on the Brownian dynamics simulations for a microliter octane droplet evaporating on glass. The velocities used in this model were found using the lubrication approximation and neglecting convective heat transfer inside the droplet [32]. Using the commercial CFD package CFD-ACE+, Dyreby et al [38] simulated the particle transport in the early stages of the evaporation of a nanoliter drop. This model assumed that the droplet shape is a spherical cap with pinned contact line throughout the evaporation, and that heat transfer and Marangoni convection can be neglected. Also, Dietzel and Poulikakos [39] simulated the coagulation of particles in a drop heated by a laser beam, a phenomenon that takes place so fast that the mass loss through evaporation is negligible. Finally, in 2008, Widjaja and Harris [40] simulated the transport of particles inside an evaporating droplet by solving a continuum advection-diffusion equation.
Due to the complex, coupled physics involved during colloidal drop evaporation, most theoretical and numerical models reported so far are based on assumptions such as fluid flow with negligible inertia [24,30,31,41], small wetting angle [19,24], spherical cap shape of the free surface [19,24,30,31,34,38,41,42], pinned wetting line throughout the evaporation [24], [30]- [35], [40,41,43], negligible heat transfer between the drop and the substrate [19,23,24,30,40,43] and negligible Marangoni convection [36,40,43]. Also, to the best of our knowledge, no published numerical study has considered the receding of the wetting line during colloidal drop evaporation, or the interaction of the free surface and the growing peripheral deposit, which can involve depinning, i.e. the separation of the drop from the deposit. Thus, a numerical modeling involving all these effects would provide a tool to predict deposit growth and final pattern shapes.
Our paper describes a numerical model that simulates the evaporation of a colloidal droplet with full consideration of the Navier-Stokes equations, convection and conduction heat transfer equations, Marangoni convection and receding of the wetting line. A detailed treatment is proposed for the interaction of the free surface with the peripheral deposit and eventual depinning. The diffusion of vapor in the atmosphere is solved numerically, providing an exact boundary condition for the evaporative flux at the droplet-air interface. The particle concentration is tracked by solving a continuum advection-diffusion equation. The model is validated against published data and used to simulate the formation of deposits during the evaporation of a nanoliter colloidal drop. Finally, the formation of different deposit patterns obtained in our laboratory is explained by our simulations.

Numerical model
The mathematical model developed in this study is based on a finite-element code for droplet impact and heat transfer developed by Poulikakos and co-workers [44]- [46] and Bhardwaj et al [47]. This model has been extensively validated for studies involving impact and heat transfer of molten metal [44,48] and water droplets [49]. It is extended here for the case of an evaporating colloidal droplet. The flow inside the droplet is assumed to be laminar and axisymmetric. All equations are expressed in a Lagrangian framework that provides accurate modeling of free surface deformations and the associated Laplace stresses [50]. Our extension of this model for the evaporation of colloidal drops involves the modeling of evaporative flux at the drop free surface, the modeling of thermocapillary stresses and the Marangoni flow, the modeling of the convection and advection of particles inside the drop and the modeling of the wetting line motion in the presence of agglomerating particles. We also propose a criterion to determine if separation of the free surface from the deposit (i.e. depinning) occurs during ring formation. Finally, we have solved the challenge to handle multiple timescales, which range from nanoseconds for capillary waves to several seconds for an entire evaporation process at ambient temperature. The various components of the modeling are presented hereafter.

Fluid dynamics
The radial and axial components of the momentum equation (equation (1)) are considered along with the continuity equation. An artificial compressibility method is employed to transform the continuity equation into a pressure evolution equation (equation (2)), as in [51,52], and a small Mach number M = 0.001 is used. Only the dimensionless form is given below, since the full derivation is available in [44].
In the above equations, M, Fr, P, τ and V are the respective Mach number, Froude number, the dimensionless pressure, dimensionless time and velocity vector, as per definitions in the nomenclature. The definitions of the dimensionless numbers are based on r max,0 = initial wetted radius of the drop and v 0 = 1.0 m s −1 as the respective reference length and velocity. The unit vector n = n R + n Z is normal to any boundary considered; n r and n z are the radial and axial components of the vector n, respectively, and T is the stress tensor: The dimensionless stress tensor terms are: No convection and radiation losses (equation (12) Balance of pressure, surface tension forces and viscous stresses (equation (7)) Hydrodynamic and thermodynamic vapor-liquid jump conditions (equation (18) and (19) Axisymmetry No slip (equation (6) Perfect thermal contact

Constant temperature
No convection and radiation losses (equation (12) No penetration for vapor concentration where Re is the Reynolds number defined in the nomenclature. It is worth mentioning that due to the Lagrangian formulation the mesh moves with the fluid velocity and advective terms do not appear in the momentum conservation equations, although they are physically considered.
The boundary conditions for the fluid dynamics are described in figure 2. Symmetry along the z-axis implies: The no-slip and no-penetration boundary conditions are applied at the z = 0 plane, except at the wetting line: Boundary conditions at the wetting line (R CL , 0) are discussed in section 2.5 because they also depend on the evaporation model. At the free surface, a stress balance is considered including forces due to surface tension, viscous stresses and thermocapillary stresses [39,53,54]: where We is the Weber number defined in the nomenclature. In the above equation,H is the dimensionless free surface curvature. The temperature dependence of the surface tension γ is assumed to be linear as γ = γ 0 + (∂γ /∂ T )(T − T 1,0 ), where subscript '1, 0' represents the initial drop temperature. The Weber number is modified to take into account the variation of surface tension with temperature: The initial conditions are that the initial shape of the drop is a sessile spherical cap with a slight overpressure due to Laplace stresses [55]: where φ 0 is the initial wetting angle of the drop.

Heat transfer
The energy equation is solved in both the droplet and the substrate, according to the formulation in [44,46,47]. Thermal convection and radiation heat transfer from the droplet free surface are neglected. A perfect thermal contact between the drop and the substrate is assumed. The dimensionless energy equations for the drop (i = 1) and the substrate (i = 2) are given by where C i is the dimensionless heat capacity and Pr is the Prandtl number. The symbol θ i is the dimensionless temperature defined as where T 1,0 and T 2,0 are the initial dimensional temperature of the drop and the substrate, respectively. In the isothermal case, a reference value of 0 • C replaces T 2,0 in the above equation.
Note that this heat transfer model accounts for both conduction and convection in the drop because of the Lagrangian approach. The boundary conditions for the heat transfer are also given in figure 2. At the substrate horizontal surface exposed to air, we neglect convection. Also, by symmetry the gradient of the temperature perpendicular to the R = 0 axis is zero. Both conditions are expressed as Along the side and bottom of the substrate (shown as BC and CD, respectively, in figure 2), a constant temperature boundary condition is applied: θ 2 (R, Z , τ ) = 1. The substrate is a square of 4 × 4 non-dimensional size, which corresponds to a dimensional value of about 1 mm × 1 mm. The initial conditions for the droplet and substrate are

Evaporation model
Since the evaporation model is part of the novel material presented in this paper, it is fully derived below. The evaporative mass flux j at the free surface of the evaporating drop is obtained by solving the diffusion equation for the vapor concentration c (kg m −3 ) [30,31,41]: where D AB is the vapor-phase diffusivity of the droplet fluid in air (for water drop, D AB = 2.6e − 5 m 2 s −1 at 25 • C [56]). Importantly, the time required for the vapor concentration to adjust to changes of the droplet shape and surface temperature is of the order of r 2 max,0 /D AB [30], which is three orders of magnitude smaller than the total evaporation time of the nanoliter drops considered in this study. Therefore, the vapor concentration evolves in a quasi-steady manner with respect to the drop evaporation and we neglect the first term of the above equation. The boundary conditions for equation (14) are also described in figure 2: where r max and z max are the wetted radius and the height of the drop, respectively. In the above equation, c vap is the saturated concentration (kg m −3 ) of the vapor and c ∞ = H c vap is the concentration in the far-field corresponding to the ambient relative humidity H. The saturated concentrations are fitted from the data in [57,58]: where T is the temperature in • C. We found sufficient to consider the far field (r = ∞, z = ∞) at r = 20r max,0 , z = 20r max,0 (figure 2). The evaporative mass flux j at the free surface can be expressed as in [31,41]: At the free surface of the drop, the hydrodynamic and thermodynamic vapour-liquid jump conditions [34,53] are applied: where v is the velocity of the liquid at the free surface and v int denotes the velocity of the free surface. The symbol n is the outward normal unit vector at the liquid-air interface and ρ is the density of the drop liquid. Also, X int is the volume concentration of the particles at the free surface and L is the latent heat of evaporation of the liquid (J kg −1 ).

Particle transport
The governing equation for the particles transport is given by [1,6,59]: where X is the concentration of the particles (kg of particles (kg of solution) −1 ) and D PL is the diffusion coefficient of the particles in the drop liquid (m 2 s −1 ). The values of the diffusion coefficient for 100 nm particles in water, 1 µm particles in water and 1 µm particles in isopropanol are calculated using the Stokes-Einstein equation [60] as 4e − 12, 4e − 13 and 2e − 13 m 2 s −1 , respectively. Boundary conditions are given as in figure 2: at r = 0, ∂ X/∂r = 0 (axisymmetry).
at z = 0, ∂ X/∂z = 0 (no penetration of the particles into the substrate).
We neglect any attraction force between the substrate and the colloidal particles because in all the systems considered here, the substrate and colloidal particles repel each other, having zeta potentials of the same polarity. We also assume that the particles do not affect the viscosity of the droplet liquid. The respective non-dimensional expressions for the hydrodynamic and thermodynamic vapour-liquid jump conditions (equations (18) and (19)) are: j · n = (V − V int ) · n and J Re Pr J a = −K 1 ∇θ · n, where J = j/ρ l v 0 and the Jacob number In the isothermal case, a reference value of 0 • C replaces T 2,0 in the denominator. The non-dimensional particle transport equation can be formulated as

Boundary conditions at the wetting line and criterion for depinning
As shown in figure 3 in our model, the drop is sessile and the wetting line is allowed to recede only if the following two conditions are simultaneously verified: • the wetting angle is smaller than the receding wetting angle (φ rec ), and • the concentration of particles at the interface (X CL ) is lower than the maximum packing concentration (0.7 for spherical particles [61]).
Otherwise, the wetting line does not move: it is pinned. In the pinned case, the Laplace stresses (equation (7)) are extrapolated to the wetting line (natural boundary condition) and the wetting line velocity v · n = v int · n + j · n/ρ(1 − X int ), is an essential boundary condition obtained from equation (18), with v int = 0 at the pinned point.
In the case where the wetting line is receding, the Laplace stresses (equation (7)) are extrapolated to the wetting line (natural boundary condition) and the wetting line velocity v is obtained as an essential boundary condition from equation (18), v · n = v int · n + j · n/ρ(1 − X int ), where v int is calculated using the fact that the wetting angle is constant (φ = φ rec ).  Receding starts when the wetting angle reaches a given receding value φ rec and then proceeds at constant wetting angle. Figure 3 explains the depinning or receding of the wetting line during the evaporation of a pure liquid drop. According to the Young-Dupre equation [29,55], depinning starts when the wetting angle of the drop φ becomes equals to φ rec (critical value of the receding angle, usually determined experimentally). As shown in figure 3, during the early stages of the drop evaporation, the wetted radius remains constant and the wetting angle φ decreases until φ rec is reached at t = t 2 . After that, depinning occurs and the wetted radius decreases while the wetting angle remains constant at its receding value (φ = φ rec ), until complete evaporation of the drop.
In the case of a colloidal drop, a particle concentration X corresponding to the maximum packing (X = 0.7) might be reached at the wetting line so that a ring starts growing as in figure 4. As the wetting line becomes cluttered by the assembling particles, we propose in figure 4 to describe the interaction of the free surface of the drops with the growing deposit as follows. First, we consider that once a closely packed ring starts forming, the wetting line shifts from its natural location (along the edge of the drop and on the substrate, at point C in figure 4) to a different location: the intersection of the ring and the free surface of the drop (point (I) in figure 4). In figure 4, the boundary of the growing ring is determined by the isoconcentration front of 0.7 shown as a thick red line inside the drop in figure 4. To propose a mechanism for depinning of the free surface from the growing particle ring, we will use the receding wetting angle φ rec (left column of figure 4) of the pure liquid drop on a substrate made of closely packed colloidal particles. This angle φ rec can be obtained experimentally and is needed to explain the interaction of the colloidal drop free surface with the growing ring (right column of figure 4), during evaporation on a smooth solid substrate. The top frames in figure 4(a) show how depinning occurs at point (I). At t = t 1 , the wetting angle φ at point (I) is greater than φ rec . As evaporation proceeds, the wetting angle φ decreases until reaching at t = t 2 the depinning value φ = φ rec . At that instant, depinning occurs and the depinning line slides along the 0.7 concentration front with φ = φ rec . The bottom frames in figure 4(b) show that a smaller wetting angle (φ rec ) for the (pure liquid)-air-particles system leads to a situation where The liquid-air-particle receding angle and its role in the colloidal ring formation process:  the colloidal ring forms without depinning because the wetting angle φ at point (I) is always larger than φ rec . Obviously, our model relies on measuring the wetting angle φ rec for the pure liquid drop on a substrate made of closely packed colloidal particle material, the value of which might be influenced by the surface roughness, of the order of the particle diameter. Note that at the interface between the drop liquid and the ring, a horizontal velocity boundary condition is applied with a magnitude matched to the evaporation rate at the surface of the ring. Also, mesh elements of which all the nodes have reached 0.7 are removed from the fluid dynamics calculation, according to the scheme defined in [62].

Numerical scheme
The mathematical model presented in the previous sections is solved using the Galerkin finiteelement method [63] with a fluid dynamics scheme from Bach and Hassager [64], as described in [44]. That scheme is, however, specifically suited for transient drop impact problems, with a time step constrained to the smallest timescale of the problem, which happens to be the oscillation of parasitic capillary waves on the droplet free surface. This timescale of x/a is about 15 ns, with x the grid size (∼2 µm) and a is the speed of sound. However, the entire evaporation process takes about 2 min for a 20 nl water drop at ambient temperature: this would   (9) and (14)- (18)) is obtained using a short time step ( τ short , order of 15 ns) with the Bach and Hassager scheme [64], and we account for evaporation with equation (18) to determine the shrinking of the free surface. Secondly, we perform a quasi-steady Lagrangian integration with a long time-step, about two orders of magnitude larger ( τ long , order of 1 µs), using the converged velocities from the previous, short time step. The computation proceeds by alternating between the two time steps. This two-step approach is justified because the fluid flow and the evaporation flux can be considered steady for the duration of the long time step because τ long τ evaporation . To compensate for the artificial loss of particles due to the shrinking of the free surface, a source term S = X jδ A/ρδV is added to the diffusion equation (equation (20)) for the elements touching the free surface of the drop, where X is the particles concentration, ρ is the density of the drop liquid, j is the evaporative mass flux and δ A and δV are the free surface area and volume of the element, respectively. Boundary conditions at the wetting line and on the ring are applied as described in section 2.5. The typical computational domain is shown in figure 2 and involves a spatial discretization of 650 nodes in the drop of linear, triangular elements, refined near the wetting line, to ensure grid and time step independency. The influence of grid size and time step have been extensively checked for heat transfer and fluid dynamics case in [44,65] for the evaporation of pure liquids. The typical computational time was about two days on an Intel dual core processor.

Thermophysical properties and dimensionless numbers
The thermophysical properties used in the simulations for water, isopropanol, glass and polydimethylsiloxane (PDMS) are given in table 1, with values taken from [56,57,66,67]. The temperature dependence of the dynamic viscosity µ, vapor-phase diffusivity of the droplet fluid in air D AB and latent heat L is taken into account. All quantities are defined in SI units, as shown in the nomenclature. For this purpose, the following curves are fitted on the data in [56,58].
For isopropanol [56,58] where the temperature T is in • C.

Code validation
While several aspects of our numerical modeling have been extensively validated, such as the fluid dynamics [68,69] and associated heat transfer [44], we evaluate below the ability of the present modeling to simulate the evaporation of colloidal drops. First, we apply our model to the evaporation of pure liquids and compare the results with published data. Then, we test the conservation of mass and energy during colloidal drop evaporation. Finally, we compare the simulated transient drop volume and accumulated mass in the ring with published analytical models of evaporating colloidal drops.  figure 5. Thermophysical properties such as ∂γ /∂ T = −1.19e − 4 Nm −1 K −1 are taken from [75]. Figure 5(a) compares the numerical and experimental transient evolution of the liquid volume and drop radius. The agreement between experiments and simulations is excellent. Also, the nonlinear evolution of the wetted radius during the receding phase is compared with a simple analytical model available in the literature [76]- [82], which can be described as follows: during a receding phase with constant wetting angle φ rec , the wetted radius for the spherical cap evolves as r max ∝ (t f − t) a , where t f is the total evaporation time. The exponent a depends weakly on the nature of the fluid, being slightly less than 0.5 for alkanes [76,80,81] but close to 0.6 for water [23,24,76].
where E 1,0 and E 2,0 are the initial thermal energies of the drop and the substrate, respectively, and E 1 and E 2 are the thermal energies of the drop and the substrate, each expressed as E i (t) = The thermal loss due to phase change at the evaporating free surface Figure 6 plots the drop volume and the relative value of the thermal energy (E 1 + E 2 + E loss ) × 100/(E 1,0 + E 2,0 ) as a function of time. Our results show that thermal energy is conserved within 1%. Also, the mass conservation is verified by plotting the mass of the colloidal particles at any time with respect to the initial particle loading (1% v/v concentration of 100 nm particles). We observe that the mass of the particles is also conserved within 1% and conclude that the modeling conserves mass and energy very well. figure 7, the temporal evolution of the particle mass in the ring obtained through the numerical model is compared with the analytical expression by Deegan et al [24], an expression valid for thin drops evaporating at ambient temperature:

Evaporation of a pinned colloidal drop. In
In equation (22), m 0 is the initial mass of particles in the drop, t f is the total evaporation time and λ = (π − 2φ 0 )/(2π − 2φ 0 ) a function of the initial contact angle. We consider two drops, 20 nl and 2 nl water drops evaporating on glass with a low value for the initial wetting angle (12 • ) and a pinned contact line, because Deegan's analysis is based on lubrication theory and is valid only for thin drops with pinned contact line. The drop and the substrate are at ambient temperature (25 • C). The comparisons of the evolution of the non-dimensional mass of the ring given by analytical theory and by the simulation are shown in figure 7. The trend of the ring mass given by the numerical results and the analytical theory are in good agreement, although the numerical growth of the ring mass for the 2 nl drop appears delayed. Figure 7 also plots the temporal evolution of the drop volume, which can be found from an analytical theory by Popov [41]: In the above equation [41,83], the drying time of the drop is proportional to the square of the value of the initial wetted radius: In equation (24), ρ l is the density of liquid (kg m −3 ), r max,0 is the initial wetted radius of the drop, D AB is the vapor-phase diffusivity of the droplet fluid in air (m 2 s −1 ), c int and c ∞ are the vapor concentrations (kg m −3 ), respectively, near the interface and in the ambient air. Figure 7 shows good agreement between our simulations and the analytical formulae above for both the mass in the ring and the drop volume evolution. However, in the later stages of the evaporation, the numerical volume decreases slightly faster than predicted analytically. A possible explanation is that our simulations were made with 5% volume of particles, whereas the correlations [24,41,83] do not specifically consider the concentration of particles. Also, results in figure 7 only show the first two-thirds of the evaporation process, corresponding to a total computation time of about two days.

Experimental details
Nanoliter water and isopropanol drops were spotted on glass and PDMS substrates, using a 375 µm diameter stainless steel pin, as routinely done to prepare bioarrays [84]. Using the same standard pin, deposited water drops were smaller (3-5 nl) than isopropanol drops (30-40 nl) probably because of higher wettability of steel by the solvent [84]. Evaporating droplets were visualized from the side using a digital camera (Pixelink, PLA 741, 1.3 megapixel up to 100 frames s −1 ) and an Optem long-distance zoom objective. Typical time and spatial resolution were 4 frames s −1 and 1 µm per pixel, respectively. We also used an Olympus IX-71 inverted microscope to qualitatively assess the structure of the flow inside the drop: the motion of fluorescent particles (1 µm or 100 nm polystyrene particles from Duke scientific Inc., CA) was recorded at frame rates of 25 frames s −1 , and resulted in figures such as figure 16 and the associated movie (available from stacks.iop.org/NJP/11/075020/mmedia). After evaporation, profiles of the deposits were measured using an atomic force microscope (Digital Instruments Dimension 3000 with Si 3 N 4 tip) in contact mode, as in figure 13. While this method worked well for deposits on glass, measurements on PDMS substrates failed: the AFM tips broke possibly because of the sticky PDMS. Therefore, the profile of deposits on PDMS was measured by the following imprinting method [85]: firstly, the PDMS substrate (Dow corning Inc.) was coated with a thick PDMS slab having a different ratio of base and curing agent, this is to facilitate subsequent removal. After curing, the thick PDMS slab was removed and washed to remove particles from the ring imprint. Then, epoxy (Loctite Inc., quick set no. 81502) was poured and cured on the PDMS slab. The positive imprint of the deposit was finally measured on the epoxy with a surface profilometer (Dektak3), as in figure 18.

Results and discussions
First, we describe two purely numerical results related to the Marangoni flow during evaporation. The first result relates the direction of the Marangoni flow loop to the ratio of substrate and droplet thermal conductivities, as per an analytical criterion (section 4.1). The second numerical result describes the effect of the intensity of Marangoni convection on the deposit shape (section 4.2). Then, we compare numerical and experimental results for two specific colloidal evaporating cases that form different patterns: Simulations were performed at ambient temperatures using the same parameters as in the experiments, with thermophysical properties and simulation parameters as in tables 1 and 2.

Direction of the Marangoni loop
Using asymptotic analysis and experiments, Ristenpart et al [27] showed that the direction of the Marangoni loop is controlled by the ratio k 2 /k 1 of substrate versus liquid thermal conductivities Case as shown in figure 8 (right column). An interesting question is therefore to determine if our numerical loops also obey the criterion in [27], and results are shown in figure 8. The idea behind the criterion in [27] is as follows.  figure 8), the heat supplied to the wetting line comes preferentially from the substrate resulting in an inversion of the thermal gradient along the drop free surface so that the Marangoni loop becomes counterclockwise. In figure 8 (left), we consider a system of 38 nl isopropanol drop on PDMS substrate with an initial wetted radius and wetting angle of 437 µm and 32 • , respectively. For this system, Ristenpart et al [27] found the critical ratio for k 2 /k 1 (when the Marangoni loop switches direction) to be 1.57. Isotherms and streamlines in figure 8 (left) show that this switching is successfully reproduced by our code, as shown by the different rotation direction for k 2 /k 1 of 1.77 and 0.6, respectively. Note that a similar reasoning can explain the counterclockwise nature of the flow in the toluene drop reported in the section 2.8.1, for which k 2 /k 1 = 2.05 is larger than a critical value of about 1.6 (see figure 5(b), right). Interestingly the numerical simulations predict a more complex flow pattern for k 2 /k 1 = 1.2 with two loops, which might occur to the proximity of the transition region. In figure 8 (left), the thermal conductivity of the substrate was the only parameter varied, all other drop and substrate properties being kept constant. about 16 s, the wetting line pins and starts forming a ring visible as an inset in figure 9 (left column). Figure 9 (right column), however, does not exhibit any ring. This numerical result shows that by reducing the intensity of Marangoni convection inside the drop, it is possible to change the final deposit pattern. Figure 10 and the associated movie (available from stacks.iop.org/NJP/11/075020/mmedia) describe our numerical results for the evaporation of a 3.7 nl colloidal water drop on a glass substrate at ambient temperature, with an initial 1% (v/v) concentration of 100 nm polystyrene particles. Time goes from top to bottom. It shows streamlines and velocity amplitude (right) as well as isoconcentrations of particles (left, labelled X). The initial value of the drop wetted radius and wetting angle obtained from experiments are 231 µm and 22 • , respectively. As in simulations from other groups involving water [32,86], Marangoni stresses are neglected because the high surface tension of water favors the deposition of contaminants at the water/air interface, a process that inhibits Marangoni convection. The parameters used in the simulation  [19] and Hu and Larson [31]. In relation to that we plot in figure 11 the evaporation flux, obtained numerically at t = 0 s, with respect to the radial distance along the substrate. Our flux results agree extremely well with the correlation of Hu and Larson [31], with both approaches showing the strongest evaporation at the wetting line.

Evaporation of a 3.7 nl water drop on glass
In figure 10, the advected particles accumulate near the wetting line and start forming a ring as soon as the particle concentration reaches the 0.7 limit corresponding to maximum particle packing, visible in red color from t = 8 s. A magnified view of the growing polystyrene particles ring is visible at the bottom of figure 10, from t = 11.7 s. At that point, the drop has flattened into a film and we have skewed the vertical scale. Since the receding wetting angle of water on polystyrene is about 80 • [87], depinning of the drop free surface from the ring might occur according to the modeling in figure 4(a), and described in section 2.5. However, our simulation in figure 10 assumes a receding angle of water on polystyrene of 45 • to account for the roughness of the ring. This prevents depinning as in figure 4(b). As a result, depinning does not occur and the ring is formed by the drying of the central film. While it is difficult to assess how the depinning angle varies with the particle diameter and associated ring roughness, we show in figure 13 that for the colloidal system considered in figure 10 numerical simulations where depinning did not occur produced a ring shape closer to the experimental ring shape. Figure 12 compares the experimentally measured and simulated evolution of the volume, wetted radius and wetting angle during the evaporation of the same colloidal water drop on  glass. In both numerical and experimental results, the wetted radius is approximately constant during the evaporation of the drop (from t = 0 to 11.7 s). This is explained by the pinning of the wetting line, as described above. The numerical results for the evolution of the drop volume show a linear decrease from t = 0 to 11.7 s, and a similar trend can be noticed in the experiment. A linear decrease in the evolution of the wetting angle can also be noticed in numerical and experimental results. The final evaporation time of the drop is around 12 s, slightly faster in the experiment than in the simulation. Figure 13 compares the profiles of the ring measured by atomic force microscopy (AFM) with the simulation result. Three AFM measurements of deposits generated under the same conditions as in figure 10 are quite similar. They show that the ring has the profile of a dome, about 1.7 µm high and 35 µm wide. Two profiles are obtained numerically, each corresponding to a different value of the depinning angle for water on the ring material, as described in figure 4 and above. The numerical profile corresponding to the lower value of the depinning angle (φ rec = 45 • , as in simulation in figure 10) appears to match the experiments quite well, better than the profile with the larger depinning angle (φ rec = 85 • ) where depinning seems to have occurred too abruptly. While both numerical simulations match the outer slope of the ring very well, the inner slope suggests that the model with φ rec = 45 • is more appropriate, and this might be explained by the fact that the roughness due to the particles assembling in the ring reduces the value of the depinning angle of the water drop.
An interesting question regarding the numerical modeling presented in this paper is to determine how large the particles can be before the continuum modeling for particle transport becomes inadequate. Therefore, numerical and experimental results are shown in figure 14 for the same case as figure 13, except that 1 µm particles are used (10 times larger than the 100 nm particles in figure 13). AFM measurements in figure 14 show a ring with a rectangular   cross-section corresponding exactly to one monolayer of particles. While the ring width is comparable to the width obtained with 100 nm particles (figure 13), the height does not compare well with the simulation results in figure 14, for the reason postulated directly above. Figure 15 and the associated movie (available from stacks.iop.org/NJP/11/075020/mmedia) describe our numerical results for a situation where the outcome is not a ring shape but a homogeneous deposit. The system corresponds to a 38 nl isopropanol drop evaporating on a PDMS substrate at ambient temperature. The drop is initially loaded with 0.1% concentration of 1 µm polystyrene particles. Figure 15 shows streamlines and velocity amplitude (right) as well as isoconcentrations of particles (left, labelled X). The initial value of the drop wetted radius and wetting angle are obtained from our experiments as 437 µm and 32 • , respectively. Marangoni stresses are taken into account with ∂γ /∂ T = −7.89e − 5 N m −1 K −1 and the physical parameters used in the simulation are given in table 2. Since the ratio of substrate and liquid thermal conductivities k 2 /k 1 = 1.77 > 1.57, we expect the Marangoni loop to turn in counterclockwise direction, according to the criterion in section 4.1. This is verified in figure 15, for t = 3-12 s. The receding angle (φ rec ) obtained experimentally is around 27 • , slightly smaller than the initial wetting angle, so that the wetting line is initially pinned and the wetting angle gradually decreases. At t = 2 s, the wetting angle hits its receding value and the wetting line recedes with a constant receding wetting angle of 27 • , according to the mechanism explained in section 2.5. Figure 15  and streamlines superposed to velocity amplitude (right) are shown. Note that the aspect ratio is increased toward the vertical axis from t = 0 to 12 s. See also the associated movie available from stacks.iop.org/NJP/11/075020/mmedia. to 12 s. As time evolves, more and more particles are advected to the top of the drop where they accumulate (t = 12 s).

Evaporation of a 38 nl isopropanol drop on PDMS
This accumulation of particles on the top of the drop agrees well with experiments in figure 16 and the associated movie (available from stacks.iop.org/NJP/11/075020/mmedia), using the same conditions as in figure 15. We visualize the flow inside the drop using fluorescence microscopy: by carefully adjusting the height of the focal plane we were able to see the radial component of the velocities at different levels in the drop, as done in [27]. The blurred image of the off-focus bead (t = 3.8 s) indicates that the circulation is counterclockwise, in agreement with our numerical result in figures 15 and 16(b). An interesting fact we noticed from our experiments is that once the drop was spotted on the substrate, the particles moved in a chaotic manner for the first 15% of the evaporation time (3 s), until a slow and axisymmetric flow took place. This phenomenon is consistent with observations made for a 2 µl isopropanol drop on PDMS substrate by Ristenpart et al [27]. The initial chaotic flow could come from oscillations that are caused by the deposition process from the pin. However, Ristenpart et al [27] speculated that this initial chaotic flow is due to the Benard-Marangoni instability. Figure 17 compares the experimental and numerical evolution of the volume, wetted radius and wetting angle for the same 38 nl colloidal isopropanol evaporating on a PDMS substrate. Experiments show that the receding of the wetting line starts at 2.5 s. The numerical results also show that a pinned phase precedes a receding phase. Also, simulations show that the wetting angle decreases linearly during the first stage of the evaporation, whereas it is constant during the second stage of the evaporation (from t = 2 to 15 s). The final evaporation time of the drop is 15 s according to the simulations, close to the measured value of 18 s.
The profile of the pattern obtained experimentally and numerically is shown in figure 18 The large oscillations and poor reproducibility of the measurement probably stem from the indirect stamping process explained in section 3. However, the general shape, volume and dimensions of the measured deposits are in good agreement with our simulation results: they show a single central bump that differs radically from the peripheral ring seen in the water-glass case of section 4.3. Our numerical results explain this difference by the change of flow pattern due to Marangoni convection, as described in figure 15.

Conclusions
In this work, the formation of deposits during the drying of nanoliter colloidal drops on a flat substrate has been investigated numerically and experimentally. A finite-element numerical model has been developed that solves the Navier-Stokes, heat and mass transport equations in a Lagrangian framework. The diffusion of vapor in the atmosphere is resolved numerically, providing an exact boundary condition for the evaporative flux at the droplet-air interface. Laplace stresses and thermal Marangoni stresses are accounted for. The particle concentration is tracked by solving a continuum advection-diffusion equation. For the first time, the receding of the wetting line and the interaction of the free surface of the drop with the growing deposit have been modeled using wetting angle criteria. Evaporation times and flow field calculated numerically are in very good agreement with published experimental and theoretical results, as well as with in-house experiments. For instance, numerical comparisons with experiments involving water and isopropanol drops loaded with polystyrene microspheres show the importance of Marangoni convection in controlling the internal flow and final deposit pattern, either a ring-like pattern or a single central bump. Numerical results are found to be in very good agreement with the measured deposit shapes and sizes.