A model for pattern deposition from an evaporating solution subject to contact angle hysteresis and finite solubility

We propose a model for the pattern deposition of the solute from an evaporating drop of a dilute solution on a horizontal substrate. In the model we take into account the three-phase contact angle hysteresis and the deposition of the solute whenever its concentration exceeds the solubility limit. The evaporating drop is governed by a film equation. We show that unless for a very small three-phase contact angle or a very rapid evaporation rate the film adopts a quasi-steady geometry, satisfying the Young–Laplace equation to leading order. The concentration profile is assumed to satisfy an advection diffusion equation subject to the standard Fick’s law for the diffusive flux. We further use an integral boundary condition to describe the dynamics of the concentration in the vicinity of the three-phase contact line; we replace an exact geometric description of the vicinity of the contact line, which is usually assumed such that mathematical singularities are avoided, with general insights about the concentration and its flux. We use our model to explore the relationships between a variety of deposition patterns and the governing parameters, show that the model repeats previous findings, and suggest further insights.


Introduction
The deposition of suspended colloidal particles and solute mass from an evaporating drop or film is a defining mechanism in coating, painting, printing, [1][2][3][4] production of low resolution electronic devices, 5 and in many scientific applications that are associated with the patterning of surfaces. 2In addition, the evaporation of particulate liquids and solutions has been studied for several decades for treating a variety of difficulties in industrial processes and developing new manufacturing procedures.[8][9][10][11][12][13][14][15][16][17] Deegan et al., 7,18,19 first pointed out that the evaporation of a drop invokes a flow field that transports suspended solid particles to the edge of the drop.Together with the pinning of the threephase contact line this gives a ring-like deposit formation -the so-called ''coffee-ring'' effect.The phenomenon was theoretically studied by Popov,8 who proposed a quasi-steady state equation for the shape of the drop using the lubrication approximation, and by Hu and Larson 20 who conducted a detailed flow simulation taking into account various mechanisms associated with evaporation, such as the temperature gradient, the vapour pressure distribution, and the Marangoni effect.The preliminary theoretical models given by Popov 8 and Hu and Larson 20 provide good approximations to the deposition process at an early stage of the evaporation.
An extension to Deegan's model was made by Fischer 21 who avoided assuming the shape of the evaporating drop is a spherical cap.Instead, Fischer proposed a dynamic shape with boundary conditions that ensure a pinned contact line.Fischer examined how different non-singular expressions for the evaporative flux affect the profile of the deposits, solving the advection equation for the concentration profile.Fisher concluded that the evaporative flux, assumed constant across the drop or enhanced near the contact line, supports an outward flow in the liquid that causes the ''coffee-ring'' effect.Otherwise when the evaporative flux primarily occurs at the center of the drop no ring-shaped deposit occurs since the solute or suspended particles are convected towards the center of the drop.
Furthermore, van Dam and Kuerten 22 proposed an extension to Fischer's study by calculating the curvature of the drop shape in order to account for large contact angles.They considered a concentration-dependent viscosity of the liquid inside the drop and incorporated the linear Fick's diffusional term into the advection-diffusion equation that governs the concentration profile.
One more extension of Fisher's and van Dam and Kuerten's models was proposed by Siregar et al., 23 taking into account the convection, diffusion, and adsorption of the solute onto the substrate.Siregar et al. 23 calculated the dynamic shape of the drop and the evolution of the deposition of the solute for both pinned and unpinned three-phase contact lines.While for a pinned contact line Siregar and coworkers observed a regular ''coffee-ring'', thus confirming previous results, for an unpinned contact line, they observed a new, mountain-like, deposit shape.These two different deposition morphologies led Siregar and coworkers to conclude that the contact line condition influences the shape of the solute deposition.To avoid singularity near the contact line they incorporated in their model a microscopic precursor film subject to a conjoining/disjoining pressure of molecular origin, based on the Frumkin-Deryugin model. 24riefly, a thin precursor film remained behind the contact line, thus helping to avoid singular corner geometry of the drop close to the contact line.
Another extension of Popov's model 8 was suggested by Ozawa et al. 2 Ozawa and coworkers solved the advection equation, assuming the evaporative flux to be spatially constant and the shape of the drop to be a quasi-steady spherical cap.An extension to Ozawa's model was proposed by Okuzono et al. 14 who treated the quasi-steady state drop in the pinned contact line case under the following assumptions: the evaporation rate remains constant in the solution phase which is 'far' from the contact line and vanishes in the gel phase, appearing in the proximity of the contact line whenever the concentration exceeds some threshold value.Furthermore, they assumed that the gel phase renders the concentration of the solute particles constant and the radial velocity vanishes.
Another model for an unpinned contact line was treated by Freed-Brown. 25Freed-Brown found the appearance of mountainlike deposits under the assumption of the quasi-steady state free surface of the drop, where most of the deposit appears near the center of the evaporated drop in the absence of repeating patterns.This study confirmed the results of Siregar et al. 23 Introducing the stick-slip motion of the contact line, Adachi et al. 26 observed the formation of circular stripe patterns, or rings from an evaporating drop.In their model, they introduced a friction force in the vicinity of the contact line due to the flow of particles towards the boundary of the drop, defining an effective mass density at the contact line and employing experimental values for the evaporation rate.They found that the mechanism of stripe formation originated in a competition between the capillary and the friction forces near the contact line.
On a different note, using the lubrication approximation, Warner et al., 6 studied the dewetting dynamics of ultrathin films in the presence of evaporation.They derived and solved two coupled evolution equations for the geometry of an evaporating film and for the corresponding concentration profile where the viscosity depends on local particle concentration.Their results indicate that starting from a uniform solute concentration the solute particles are aligned to distinct deposits, where the spacing between the deposits depends on the parameters of the evaporating system.
Avoiding empirical conditions on the contact line and the contact angle, Thiele et al. 9,10,27 have derived and solved two coupled evolution equations for the geometry of an evaporating film and for the corresponding concentration profile, assuming the presence of a disjoining/conjoining pressure of molecular origin and a precursor film.Thiele and coworkers employed the standard Fick's law, where the diffusion is dependent on the excess Gibbs free energy of the mixed phases in the evaporating film, and a power relationship between the viscosity and the concentration.They showed that each of the non-linear contributions to diffusion and viscosity is independently sufficient to supporting pattern deposition.
An additional model for the stick-slip motion of the contact line in the case of a moving substrate was suggested by Doumenc et al. 28,29 Doumenc and coworkers solved two coupled evolution equations for the geometry of an evaporating film above a moving substrate and for the corresponding concentration profile.The viscosity was assumed to exhibit a strong relationship with the volume fraction of the solute and was prescribed as an empirical function of solute concentration.They assumed that outside the meniscus the system is frozen due to the very high viscosity of the solution, which is analogous to the assumption of the precursor film.Doumenc and coworkers obtained periodic patterns whose amplitude was found to depend on the velocity of the moving substrate.Furthermore, Doumenc et al. 30,31 have found an empirical law that relates the pinning force and the deposits' geometry (in the case of a moving substrate), and concluded that this force controls the wavelength of the periodic patterns.
A piece-wise model that accounted for the stick-slip motion of the contact line was proposed by Kaya et al. 32 Kaya and coworkers solved the advection-diffusion equation, where avoiding singularity the dynamic assessment of the concentration was performed in the immediate vicinity of the contact line.The shape of the drop was assumed to be a spherical cap and the contact line was assumed to unpin when the maximal concentration of the suspended particles exceeds an empirical threshold value.Contact line pinning was assumed to occur at an empirical equi-spaced distance between the concentric rings.
Furthermore, studying the stick-slip motion of the meniscus in the case of a Petri dish geometry, Chen et al. 33 proposed a piece-wise model that takes into account the antagonistic balance between the pinning and de-pinning mechanisms, assuming a constant evaporation rate.Chen and coworkers related pinning/de-pinning mechanisms, the receding contact angle, and the roughness of the substrate, where they based their analysis on the relationship between the receding contact angle and the roughness of the solid substrate obtained by de Gennes. 34he simulations of Chen and coworkers have qualitatively yielded a quasi-logarithmic spatial relation between the deposits.
In this paper, we present a piece-wise model, reminiscent of the models by Kaya et al. 32 and Chen et al., 33 for the pattern deposition of the solute mass from an evaporating drop whose contact line undergoes a stick-slip motion as a function of the receding contact angle threshold on the bare substrate and on the deposit.We assume a diffusion limited evaporation 7,8,18,19 and further incorporate the standard Fick's diffusion law 27 in the advection-diffusion equation for the diffusive flux in the solution.Moreover, we incorporate in our model the solubility limit of the solution.The solubility limit determines the rate of solute deposition to occur when the concentration exceeds the solubility limit, and the flux of concentration in the liquid, associated with dissolved solute mass.We treat the vicinity of the contact line using an integral boundary condition to account for the local dynamics of the concentration.
We use our model to study the pattern deposition for pinned and un-pinned contact lines and for an interchanging stick-slip motion of the contact line.We characterise the patterns by performing a parametric study, confirming previous experimental results 7,16,18,19,33,35 and giving further insights.In the following we derive our model in Section 2, discuss the method of solution in Section 3, summarize our findings in Section 4, and conclude in Section 5.

Model
In this section we derive a model for the pattern deposition of solute mass from an evaporating drop of a dilute solution.We assume that the spatiotemporal dynamics of the drop's shape and that of the concentration may be described using the principles of lubrication theory.We ignore gravitational contributions and assume that the density r, the viscosity m, and the surface tension g of the drop are constant to leading order.
The analysis of the corresponding dynamic film equation, which we present in Appendix A1, suggests that the ratio of the capillary velocity to the rate of evaporation, given by A gy 0 4 /mJ 0 , determines the spatiotemporal geometry of the drop, where y 0 is the receding contact angle of the solution with the substrate and J 0 is the characteristic rate of evaporation per unit area.We take the characteristic rate of evaporation to be J 0 = D v c v (1 À H)/R 0 r, assuming that the rate of evaporation of the drop is limited by the rate of diffusion of the vapour in the air above the drop; 8,19,20,23,32 the rate of diffusion is determined by the vapour coefficient of diffusion in air D v , the saturated vapour concentration c v , the ambient vapour concentration H, and the initial radius of the drop R 0 .The non-dimensional parameter A is usually large, unless for a very small contact angle or a very large rate of evaporation.For example, an evaporating drop of water under ambient conditions gives A E 10 5 -10 6 .Thus, to leading order the shape of the drop is governed by capillary stress.The drop then takes a quasi-steady hemispherical shape, changing in time due to evaporation.An alternative justification for the hemispherical equilibrium drop's shape was derived by Popov et al. 36 We assume a radial symmetry for the sessile drop and use the cylindrical coordinates (r,z).We denote the time-dependent drop's projection radius on the flat substrate by R(t), where t is time.Furthermore, the free surface of the drop satisfies the three phase contact angle y(t) with the substrate, see Fig. 1.For simplicity, throughout this work, we shall use the following notation for the partial derivatives: q r q/qr, q t q/qt, etc.We start our discussion by considering two limiting cases of the stick-slip motion of the contact line, the pinned and the unpinned contact lines, which were widely considered in the literature. 2,7,8,23,25Assuming large A, the thickness of the drop h(r,t) in these two limiting cases satisfies to leading order the following equations and boundary conditions: ð0; RÞ; @ r hð0; tÞ ¼ hðR; tÞ ¼ 0; @ r hðR; tÞ ¼ Ày; where Dp is the (unknown) Laplace pressure of the drop, which may depend on time, and dV/dt is the integral volume loss of the drop due to evaporation, to be expressed explicitly in the sequel.In the pinned contact line case the position of the contact line is fixed.However, the contact angle y = y(t) decreases in time.In the unpinned contact line case the contact angle is assumed to be independent of time, taking the value of the receding contact angle.The radius of the drop R = R(t) decreases in time.Note that the equation in ( 1) is a second order ODE.In addition, Dp and either y(t) or R(t) are unknown.
In the first case y 0 = y(t = 0) is prescribed whereas in the second case R 0 = R(t = 0) is prescribed.Thus, in both cases, problem (1) results in four unknowns and four constraints.The problem is solvable analytically.We give explicit expressions for the dimensionless shape of the drop in the pinned and un-pinned contact line cases in Appendix A3.
The solution to the problem in (1) provides us with the shape of the free surface of the drop.The same problem as the problem stated in (1) was solved by Popov, 8 with R(t) and y(t) both depending on time simultaneously.Hence, Popov needed to make an additional assumption with regard to the vicinity of the contact line that related between R(t) and y(t).Furthermore, a similar problem for h(r,t) was solved by Ozawa et al., 2 in the pinned contact line case, under an additional assumption that the evaporative rate is constant, independent of time.
Introducing the equation for mass conservation in the drop 2,8,32 q t h = À(q r /r)(rhu ¯) À J(r)/r, we further extract the average of the radial velocity uðr; tÞ ð1=hÞ Ð z¼h z¼0 uðr; z; tÞdz, where, u(r,z,t) is the corresponding velocity along the same direction and J(r) is the local evaporative flux.Integrating from the axisymmetric center of the drop and rearranging the equation gives providing us with the convectional velocity of the solute in the drop to result from the evaporative flux.Furthermore, Deegan et al. 7,18 highlighted the analogy between the fields of vapour concentration and electrostatic potential for estimating the rate of evaporation.Analytical solutions to the electrostatic potential field are given by Lebedev 37 for various geometries.The local rate of evaporation calculated from a vapour concentration field, corresponding to an evaporating hemispherical drop of a small contact angle, is 8,19,20,23,32 JðrÞ ¼ 2 p where c s is the saturated vapour density and c N is the ambient air density.It then follows that since the only contribution to the mass loss from the drop is evaporation, the integral mass loss of volume is given by where V ¼ 2p Ð r¼R r¼0 rhdr is the time-dependent volume of the drop.The local average concentration in the drop in terms of density fðr; tÞ ð1=hÞ Ð z¼h z¼0 cðr; z; tÞdz is governed by the film concentration equation, 32 where the diffusive flux is governed by the standard Fick's law, 27 such that where c, D, and f B (f/r) E (f/r)ln(f/r) + (1 À f/r)ln(1 À f/r) are the local concentration of the solute in terms of density, diffusion coefficient, and the excess Gibbs free energy per unit volume for a dilute mixture (divided by rRT), respectively, where R is the ideal gas constant, T is the absolute temperature; and f B 0 is the complete derivative of f B .For a detailed derivation of ( 5), see Appendix A2.
To render the equations non-dimensional, we use the transformations r -R 0 r, h -R 0 y 0 h, yy 0 y, t -(R 0 y 0 /J 0 )t, frf, J -J 0 J, and u ¯-( J 0 /y 0 )u ¯, where R 0 is the initial radius of the drop.
The dimensionless form of the advection-diffusion equation in ( 5) is given by Thus, the concentration profile is governed by the ratio of the diffusion velocity to the rate of evaporation D Dy 0 /J 0 R 0 .
In order to solve the equation we prescribe the homogeneous initial concentration f(t = 0) = f init , the symmetry condition q r f| r=0 = 0 at the center of the drop and an integral boundary condition for the concentration in the vicinity of the contact line near r = R.We derive the integral boundary condition in the next subsection.

Integral boundary condition
The singular corner geometry of the drop close to the contact line gives rise to a diverging evaporative flux and a diverging velocity field. 7,8,19,20Different approaches for eliminating the singularity in eqn (6) in the proximity of the contact line were discussed in the literature. 38The most common methods are: smoothing the evaporative flux singularity mathematically, 21 finding analytical solutions under simplifying assumptions, 2,14,25 cutting the proximity of the contact line, 8,32 and adding a thin precursor film which regularises the singular wedge geometry near the contact line. 9,23,39ince the mathematical singularity of the concentration, appearing near the vanishing thickness of the drop at the quasi-static contact line, is integrable we propose to treat the advection-diffusion equation in the vicinity of the contact line using an integral approach.We integrate the advection-diffusion equation in ( 6) from r = r* to R, where Dr*/R { 1 (Dr* R À r*), subject to the physical insight that the flux of concentration cannot leave the liquid and thus the drop.Only the deposited solute that left the liquid phase may leave the drop under the contact line; this is discussed in the following subsection.In order to attain an analytical solution for the concentration at r*, we further assume that changes in the mass of solute within the region r A (r*,R) are small such that the integral quantity of hf is in a quasi-steady state in this region with respect to the rest of the drop.The meaning of the last assumption is visualized in Fig. 2. The advantage of this method is that we use general insights about the concentration and its flux near the contact line rather than specific assumptions about the shape of the drop in the vicinity of the contact line.Furthermore, the results of the model appear to be not very sensitive to the magnitude of Dr* as long as 0 o Dr*/R { 1.The corresponding sensitivity analysis is given in Appendix A7.
The integral boundary condition gives the rate of change of the solute mass in time at r*: The small integral quantity of the solute mass at r* o r o R is given analytically.The condition in ( 7) is used to determine f(r*,t); all other quantities are known.Further details are given in Appendix A4.

Contact angle hysteresis and solubility limit
In our model we take into account two physical parameters that further affect the pattern deposition of the solute: the receding contact angle over the deposit y r , which differs from the receding contact angle of the drop over the bare substrate, prior to the deposition of solute particles, 16 y 0 , and the solubility limit f 0 of the solution.First, we discuss the receding contact angle y r and the mechanisms which govern the contact angle hysteresis.Then we discuss a method which allows incorporating the solubility limit in the advection-diffusion equation, so that the resulting equation constitutes a regular PDE for f.
Pattern deposition of a solute mass from an evaporating solution may invoke two different types of contact angle hysteresis mechanisms.The first is the contact angle hysteresis over the bare substrate.The second is the contact angle hysteresis over a deposited layer of solute mass.In this model we account for the two corresponding receding contact angles.One can assume that the initial contact angle of the evaporating drop, prior to the deposition of solute mass, is equal to, or greater than, the receding contact angle on the bare substrate y 0 .Motion of the contact line occurs only when the contact angle is smaller than this threshold, i.e., y o y 0 .Furthermore, once the deposition of the solute occurred in the vicinity of the contact line the chemistry of the surface at this position changes and thus the value of the receding contact angle changes to y r .
The necessity of y r in the model and the assumption that y r r y 0 may be explained as follows: if there was no y r , which reflects the receding contact angle over the deposited material, or if we had assumed that y r 4 y 0 , then we would get that y(t) o y 0 or y(t) o y r , respectively, for all t 4 t ˜, where t ˜denotes the first time for which depinning has occurred over the area of the substrate covered by deposited solute particles.So in these two cases, nothing could stop the motion of the contact line and thus no more pinned contact line stage could be achieved, meaning that no periodic patterns can be expected.
To be consistent with our previous discussion, the scaled receding contact angle over the deposit is given by the transformation y ry 0 y r .We further assume that in a scaled notation 0 r y r r 1, which also corresponds to a suggestion made by Bhardwaj et al., 16 based on experimental evidence.The two limiting cases y r = 0 and y r = 1 correspond to a pinned (once deposition occurred) and an unpinned contact line, respectively, (our calculations commence assuming that the contact angle is at its receding threshold on the bare substrate, y(t = 0) = 1).
The underlying assumption of our model is that the contact line of the drop is pinned unless the dynamic contact angle is less than the corresponding receding contact angle on the bare substrate or less than the corresponding receding contact angle on the deposited material; in which case de-pinning occurs.The details of the implementation of this approach in the numerical model are given in Section 3.
We now incorporate the solubility limit of the solution in our model.The solubility limit is a thermodynamic property of the solution, where solute particles leave the solution and deposit as an independent phase whenever the concentration f(r,t) exceeds this limit.The deposited particles do not participate in the concentration (solute mass) flux in the advection-diffusion equation.Note that our notion of the solubility limit is roughly analogous to the ''gelation concentration'' used in the model of Okuzono et al., 14 in the sense that as far as the concentration exceeds some threshold in some region (in the vicinity of the contact line), it becomes constant (equal to that threshold) in that region.
Incorporating the solubility limit implies that the flux in eqn ( 6) is determined by the component of the solute that is soluble, given by the concentration j(r,t).The actual concentration j(r,t) is limited by the solubility limit f 0 .The following modified advection-diffusion equation takes into account the total solute mass (on the left hand side of the equation), but only the mass of the soluble component of the solute when considering the convective and diffusive fluxes of the solute mass (on the right hand side): Now the 'generalised' concentration f(r,t) accounts for the overall local solute mass in both the liquid (as a solution) and outside the liquid (as a deposit).We define j(r,t) using a smooth approximation for the minimum function min{f(r,t),f 0 } that gives the minimum value between f(r,t) and f 0 , so that its value cannot exceed the solubility limit.See Fig. 2 for an illustration.Note that since j(r,t) is a smooth function of f(r,t), the equation in ( 8) constitutes a regular PDE for f(r,t).Thus, the numerical scheme which is described in the next section may converge to its solution for a given initial value of the concentration f init .See Appendix A6 for further details regarding the definition of j and Appendix A7 for the sensitivity analysis.
Modifying the boundary condition in (7) to incorporate the solubility limit (since it is derived from the advection-diffusion equation) gives The solubility limit further determines the solute mass per unit area, M(r,t), which was deposited until time t behind the contact line, and is defined as M(r,t) h(r,t*)max{f(r,t*) À f 0 ,0}, for r Z R(t), (10)   where h(r,t*) and f(r,t*) give the local thickness of the drop and solute concentration just before de-pinning occurred at t*, respectively, in order to convert the local excess concentration f(r,t) À f 0 to the mass of deposited solute that is left behind the receding contact line.See Fig. 2 for an illustration.This analysis takes place between the contact line at time t and the initial position of the contact line at time t = 0 at R(t) o r o R 0 .The function max{f(r,t) À f 0 ,0} gives the maximal value between f(r,t) À f 0 and 0, such that no deposition occurs where locally f(r,t) r f 0 .We further assume that the deposited mass of the material does not change during the de-pinning process since no net flux of volume passes through the contact line in either direction at a laboratory frame of view.Furthermore, note that the definition in (10) of solute mass per unit area is analogous to the notion of 'dry-deposit thickness', which was used to study the properties of the deposited solute patterns by Doumenc et al. 28,29 3 Numerical procedure In a similar manner to Ozawa et al., 2 and Schwartz and Eley, 24 we have used the alternating-direction-implicit technique (ADI) to step forward in time, 40 where the nonlinear prefactors in the equation are approximated at the previous time level.
The simulations subject to contact angle hysteresis are performed as follows: we assume that the initial contact angle is the receding contact angle on the bare substrate, where y(t = 0) = 1.We evaluate the initial shape of the drop, h(r,t), assuming that the contact line is pinned.See Appendix A3 for further details.At each time step we solve the equation for f(r,t) as prescribed in (8).Since y(t) decreases with time and max {0oroR} f(r,t) increases with time, the simulated drop reaches a time, which we denote by t n nÁDt, where n A N and Dt are the time-step size, for which exactly one of the conditions below is satisfied: While the contact line is pinned the maximal concentration is achieved in the proximity of the contact line. 7,18The meaning of (11a) is that the contact line should slip to the next pinned position because no deposit occurs on the substrate at this stage and the contact angle is less than the receding contact angle on the bare substrate.The meaning of (11b) is that the contact line should slip to the next pinned position because a deposit was formed in the proximity of the contact line at this stage and the contact angle of the drop is smaller than the receding contact angle on the deposited area.If one of the conditions (11a) or (11b) is satisfied, then the contact line de-pins.Otherwise the contact line is pinned.Should de-pinning occur, we determine the new position of the contact line, where the contact angle is to undertake its receding threshold on the bare substrate, subject to mass conservation.We then continue our calculations from the new position of the pinned contact line.
Let us denote by t n and R n the last time level before the de-pinning and the corresponding position of the contact line, and by R n+1 the new pinned position at time t n+1 after the de-pinning has occurred.Accordingly, we denote by y n and y n+1 the contact angles at the pinned positions R n and R n+1 , respectively.In a similar manner to Chen et al., 33 we employ the assumption that y n+1 is equal to the receding contact angle on a bare substrate, that is y n+1 = 1.Using the last constraint which was prescribed in problem (1) for the two corresponding drops' shapes (before and after the de-pinning), where the corresponding volume loss in each case is calculated according to (4), we get the following cubic equation for R n+1 : The real root of ( 12) may be easily found, thus R n+1 is determined.For further details, see Appendix A5.We may summarize our numerical procedure as follows: the number of steps in time of the algorithm is given by n.We start with an initial condition, where n = 0, R n=0 = 1, and y n=0 = 1.
(1) For n Z 1, we calculate y n y(nDt) and the corresponding parabolic free surface shape h n h(r,t n ) for the prescribed value R n R(nDt) = R nÀ1 , assuming a pinned contact line and using eqn (34) and (35) in Appendix A3.
(2) We calculate u ¯n using the formula in (2) and f n+1 using eqn (8), where f 0 = f init .We employ the no flux boundary condition for f at r = 0, and the integral boundary condition (9) at r = r*.
(3) We check if any of the conditions (11a) or (11b) is satisfied.If yes, we proceed according to step (4).Otherwise we adjust n forward by 1 and move back to step (1).
Eqn (34) gives h n+1 .Next, we adjust n forward by 1 and move to step (2) with the new y n , R n , and h n values.
We denote by t = T the end time of the simulation, where the solution has almost fully evaporated such that the drop vanishes and our solution for its leading order geometry becomes singular.Furthermore, we use the following discrete form to calculate the deposited mass M(r,T) defined in (10), where {R n } N n=0 denotes a sequence of contact line positions, starting from the initial position of the contact line R 0 and moving towards the center of the drop, and {t n } N n=0 denotes a sequence of the corresponding last time levels before the de-pinning occurs, where N is the maximal integer satisfying t N r T. Note that this definition implies that whenever the deposit M(r,T) is generated between the old pinned position R n and the new pinned position R n+1 , it is calculated as the product h(r,t n )[f(r,t n )Àf 0 ], where t n is the last time level corresponding to the pinned position R n .

Results
In this section we present our results obtained for various choices of the parameter values for the scaled diffusion coefficient D, the ratio of the receding contact angle on the deposit to the bare substrate y r , and the ratio of the solubility limit to the initial concentration b f 0 /f(t = 0).First we give representative deposition patterns.Then we employ general characterisation methods in order to compare different pattern depositions that we obtain in our parametric study.
In Fig. 3 we present the profile of the deposited solute mass per unit area on the substrate M(r,t = T) following the full evaporation of the drop.In the figure we observe the generation of detached patterns far from the center of the drop.The patterns merge into a corrugated solid film near the center.The change in the deposition state is due to the increase in concentration and capillary stress as the drop evaporates and shrinks in volume.This spatial change in the state of pattern deposition is commonly observed in experiment. 13,26,35,41Furthermore, a non-monotonic behavior of the rings' heights has been previously observed in experiment as well. 16n the pinned contact line case (which is obtained for y r = 0) we get a ''coffee-ring'', 7,18,19 which forms due to the evaporative flux invoking a flow field and thus transporting the suspended solid particles to the edge of the drop.See Fig. 4 for example.Furthermore, confirming previous findings, 23,25 in the unpinned contact line case (which is obtained for y r = 1), we observe a mountain-like deposit shape.See Fig. 5 for example.Here solute particles which are accumulated in the proximity of the contact line are gliding together with the contact line towards the center of the drop.
Next we compare between the geometries of the different pattern depositions of the solute mass M that appear under different physical and initial conditions using their relative roughness, c, the individual distances between rings l n , and the ratio of the area on which pattern deposition occurred to the initial area of the drop on the substrate P.
We measure the roughness of mass, deposited over the substrate, using the notion of 'RMS', i.e., the root mean square  operator, [42][43][44] which is one of the most common methods to define roughness.We thus employ the 'RMS' operator on the topography of the deposit.It represents the standard deviation of the variation in the mass deposited over the substrate.Furthermore, in our study, we wish to compare one deposition system to another, which may differ by the overall solute mass, using the relative roughness c of the deposited mass M(r,T), which is given by dividing the 'RMS' measure for the roughness by the spatial average of the deposited mass.The solution of the governing equations gives the 3D topography of the deposit in the form of a vector.Assuming K numerical grid points, the vector is: (2pr 1 Dr 1 m 1 , 2pr 2 Dr 2 m 2 ,. ..,2prK Dr K m K ), in which r k is the distance between the grid point k and the origin, Dr k is the distance between the grid points k and k À 1, and m k is the thickness of the deposit near the grid point k.We then divide the result by the arithmetic mean of the film thickness P K k¼1 2pr k Dr k m k .In summary, we define the roughness using Recalling that the initial radius of the drop satisfies R 0 = 1, we note that in the case of equi-spaced grid, Dr k 1/K, which allows simplifying (14), in the following manner, We visualise the meaning of the relative roughness for pinned and unpinned contact lines and for the stick-slip motion of the contact line in Fig. 4, 5, and 6(a and b), respectively.In the pinned and unpinned contact lines the relative roughness reflects the thickness of one deposited ring (in either the edge of the drop or near the drop's center, respectively), relative to the thickness of a smoothly deposited layer of the same mass, where higher values of c correspond to a higher relative thickness of the deposited ring.In a similar manner, under a stick-slip motion of the contact line, the relative roughness reflects the average thickness of the deposited rings relative to the thickness of a smoothly deposited layer of the same mass, with higher c values corresponding to thicker deposited rings.
Furthermore, for the case of a stick-slip motion of the contact line we calculate the individual distances between the N deposited rings.The distance between ring numbers n and n + 1, counted from the outermost ring, is given by l n r max n À r max n+1 , such that n = 1, 2, . .., N À 1; the coordinate r max n denotes the r-coordinate of the n local maximum value of the deposit mass M(r,T) (the maximum of ring number n).In the results we observe a logarithmic decline in the individual distances between rings l n versus the ring number n counted from the outermost ring.See Fig. 6(c and d) for example.This is the consequence of the radial geometry of this problem.This result agrees with the findings reported by Chen et al. 33 Note that the slope of the natural logarithm, ln(l n ) versus n, depends on the choice of the parametric values for D, y r , and b.Moreover, in Fig. 6 and 9(b) we demonstrate that the slope of the natural logarithm of the distance between the deposits, ln(l n ), versus n reflects the spatial density of the deposits and that the denser the deposited rings are, the larger the slope of ln(l n ) versus n becomes.In addition, P denotes the area on which the pattern deposition is apparent with respect to the initial area of the interface between the drop and the substrate.This quantity gives the scaled position of the outermost deposit (with respect to the center of the drop) and is illustrated in Fig. 6(a and b).We have performed a parametric study in three cases, stickslip motion of the contact line, pinned and unpinned contact lines, see Fig. 7-9.In all cases, the ranges for D, y r , and b were chosen such that we were able to ensure the convergence of the numerical scheme.For an example of the numerical sensitivity analysis, see Appendix A7.In Fig. 7(a and b) and 8(a), we present a graph of the roughness versus the scaled diffusion coefficient D for pinned and unpinned contact lines and for the  This difference in the trends of roughness versus D may be explained as follows: for the stick-slip motion and the pinned contact line, the contribution of the diffusion opposes the convective flux near the contact line.The convective flux tends to bring the solute particles towards the contact line while the diffusion decreases gradients in the concentration.Thus, decreasing the contribution of the diffusive flux by decreasing D supports higher 'peaks' of concentration near the contact line, which translates into a larger roughness of the deposit.
However, for the unpinned contact line the dynamics is different.While the flow in the drop convects solute particles towards the contact line, the contact line moves towards the center of the drop, such that the convective flux vanishes under ideal conditions when considering a laboratory frame of view.The concentration of the solution accumulates now due to the evaporation of the drop.The accumulation of the concentration is most pronounced near the moving contact line and as before it is dispersed towards the center of the drop by the diffusive mechanism.Decreasing D, we get a ''smaller push'' of the concentration towards the center of the drop, which results in a lower concentration peak in the center of the drop and more solute mass being deposited behind the moving contact line.This translates into a lower roughness of the deposit with decreasing D. Furthermore, as shown in Fig. 7(a and b) and 8(a), the roughness is a monotone increasing function of b in all cases.
The graph of the roughness versus the receding contact angle y r of the deposit is presented in Fig. 9(a), where again  the various trend-lines correspond to different values of b and, following our previous discussion, the receding contact angles y r = 0 and y r = 1 correspond to the pinned and the unpinned contact lines, respectively.We may conclude that the roughness is a monotone decreasing function of y r and again it is a monotone increasing function of b in all three cases, as previously.
The graph of the position P of the outermost deposit versus scaled diffusion coefficient D is shown in Fig. 8(b), where the trend-lines represent again different values of b.Fig. 8(b) predicts that increasing either D, or b, we obtain that the position P of the outermost deposited peak is closer to the center of the drop, giving a smaller deposited area.
Lastly, the upper panel of Fig. 9(b) presents the slope of ln(l n ) as a function of n versus the receding contact angle y r of the deposit and the lower panel of Fig. 9(b) presents the natural logarithm of the number of deposited rings versus y r .From these two panels we may conclude that the slope and the density of the deposited rings, which are shown to be equivalent, are mainly affected by the receding contact angle of the deposit.Increasing the receding contact angle of the deposit, we observe more peaks in our pattern, which become denser where the slope of ln(l n ) increases as well.According to Fig. 9(b) it appears that the density of the deposits is approximately independent of b, since the various trend-lines approximately coincide in both panels.

Conclusions
In this study, we have proposed a model for pattern deposition from an evaporating solution subject to the stick-slip motion of the contact line and subject to the two limiting conditions of the stick-slip motion, given by the pinned and the unpinned contact lines.We assume that solute particles are deposited whenever the concentration exceeds the solubility limit and we propose to treat the mathematical singularity of the concentration and the volume flux near the quasi-static contact line using an integral boundary condition.We obtain qualitative agreement between the geometry of the deposited mass obtained in previous experimental and theoretical studies 2,16,18,19,33,35 and our results for pinned and unpinned contact lines and for the stick-slip motion of the contact line.
In the case of the stick-slip motion of the contact line our model predicts the deposition of concentric rings.The lengths of the spacings between the rings decrease towards the center of the evaporating drop in a logarithmic manner.To further characterise the obtained patterns we extracted the roughness of the deposits (giving a measure of the geometrical deviation of the deposits from a corresponding smooth layer of deposit of the same mass), the variations of the spacings between consecutive rings (giving a measure of the spatial density of the deposited rings), and the position of the outermost deposited ring with respect to the initial radius of the sessile drop (giving a measure of the effective area for deposition), from our numerical data.
Following a parametric study we observe that for the case of a stick-slip motion of the contact line, roughness is a monotone decreasing function of the scaled diffusion coefficient D and of the ratio between the receding contact angle on the deposit and the bare substrate y r .Furthermore, the roughness is a monotone increasing function of the ratio between the solubility limit and the initial concentration b.
For the pinned contact line we observe the same trends of roughness versus D and versus b as in the case of the stick-slip motion.For the unpinned contact line, we observe a monotone increase in the roughness versus b, however we obtain an opposite trend of roughness versus D. Namely, for the unpinned contact line the roughness is a monotone increasing function of D for small values of b, whereas for high b levels the roughness appears to slightly decrease following a maximum value.
The position of the outermost solute deposition with respect to the center of the drop, which reflects the overall deposited area, is dictated by the diffusion coefficient D and the ratio b.Increasing either D, or b, we obtain that the outermost deposited peak is closer to the center of the drop, giving a smaller deposited area.
The deposition-density of the rings, which was shown to be equivalent to the slope of natural logarithm of the individual distances (between the deposited rings) ln(l n ) versus the number of the individual deposited ring, commencing from the most outer ring, n is mainly affected by the ratio of the receding contact angle on the deposit to the bare substrate y r .Increasing y r we observe more peaks of deposited solute in our pattern, which become denser.Furthermore, the density of the deposits is approximately independent of b.
This journal is © The Royal Society of Chemistry 2016 q r (rh 3 q r (q r /r)(rq r h)) = (3r/A)( J(r) À q t h) suggests that to O(A À1 ) the geometry of the film h(r,t) is quasi-static and satisfies q r (rh 3 q r (q r /r)(rq r h)) = 0, excluding the vicinity of the three phase contact line where the rate of evaporation becomes large.Integrating this equation twice and requiring that h(r,t) is not singular gives the Young-Laplace equation in Appendix A3.

A2 The advection-diffusion equation
Let us denote the fluid-air interface by h(x,y,t) over a rectangular domain (x,y) A (0,L) 2 , where t A (0,T) denotes time.We start from the basic equation for motion in 3D: q t c + rÁ(uc) = ÀrÁJ, (x,y,z) A (0,L) 2 Â (0,h(x,y,t)), (17a) where c(x,y,z,t) is the concentration in terms of density, u = (u(x,y,z,t), v(x,y,z,t), w(x,y,z,t)) is the velocity field, r (q x ,q y ,q z ), n is the unit normal to the surface z = 0 or z = h(x,y), and q n is the normal derivative to the corresponding surface.Eqn (17a) is known as the continuity equation for a mixture, and the condition in (17b) is a consequence of no penetration boundary condition, since no penetration boundary condition implies that no solute particles are allowed to penetrate neither through the substrate nor through the free surface of the liquid.If the flux, J, is defined by J Àc D RT rm, where D is the diffusion coefficient, R is the universal gas constant, T is the absolute temperature, and the chemical potential for one component in a solution is given by m m 0 + RT ln(c/r), where r is the density of the solution as defined previously, then eqn (17a) may be expressed as, where D is the Laplacian operator.
To render the equation non-dimensional, we use the transformations xl 0 x, yl 0 y, zl 0 y 0 z, hl 0 y 0 h, tt 0 t, u -(l 0 /t 0 )u, v -(l 0 /t 0 )v, w -(l 0 y 0 /t 0 )w, and crc, where l 0 is the initial length of the rectangular domain, 0 o y 0 { 1 is the receding contact angle of the bare substrate, and t 0 is the characteristic time scale.
Considering the boundary conditions which were prescribed in (17b), along the boundary z = 0, we get that q z c 0 (x,y,0,t) = 0.
From ( 21), (22), and ( 24), we conclude that at leading order @ z 2 c 0 ðx; y; z; tÞ ¼ 0; ðx; y; zÞ 2 ð0; LÞ 2 Â ð0; hðx; y; tÞÞ; @ z c 0 ðx; y; 0; tÞ ¼ 0; ðx; yÞ 2 ½0; L 2 ; @ z c 0 ðx; y; hðx; yÞ; tÞ ¼ 0; ðx; yÞ 2 ½0; L 2 : From ( 25), it follows that c 0 (x,y,z,t) is independent of z, or in other words, that there exists a function f(x,y,t), so that c 0 (x,y,z,t) = f(x,y,t), where integrating the expansion in (20)  with respect to z, 0 r z r h(x,y,t), we get at leading order that fðx; y; tÞ ¼ 1 h Ð h 0 cðx; y; z; tÞdz.Using the fact that the leading order term of c(x,y,z,t) (denoted by f(x,y,t)) is independent of the z-coordinate and multiplying eqn (19) by y 0 2 , we get at order O(y 0 2 ) that q t f + q x (uf) + q y (vf) + q z (wf) = D(q x 2 f + q y 2 f + q z 2 c 1 ), (26)   where we define D Dt 0 /l 0 2 .Integrating the equation in (26) with respect to z, 0 r z r h(x,y,t), and employing the Leibnitz rule, we get that where uðx; y; tÞ 1 h Ð h 0 uðx; y; z; tÞdz and vðx; y; tÞ 1 h Ð h 0 vðx; y; z; tÞdz are the averaged velocities in the x and y coordinates, respectively.Note that by the no-penetration boundary condition it follows that w(x,y,0,t) = 0.Moreover, from the boundary condition in (17a), we get at order O(y 0 2 ) that q z c 1 (x,y,0,t) = 0. ( As to the boundary condition at z = h(x,y,t), from (23), we get at order O(y 0 2 ) that (Àq x h,Àq y h,1)Á(q x f,q y f,q z c 1 ) = 0, at z = h(x,y,t), which implies that q z c 1 = (q x h)(q x f) + (q y h)(q y f), at z = h(x,y,t).(29)  This journal is © The Royal Society of Chemistry 2016 and R n+1 is the new position at time t n+1 after the de-pinning has occurred.From (34), we get that h Thus, during the de-pinning the drop loses the following amount of volume, where Dt t n+1 À t n .We may determine the new position of the contact line analogously to Chen et al., 33 although in ref. 33 the underlying assumption is that the volume of the drop remains unchanged during the slipping process.However, we assume that the liquid continues to evaporate during the slipping process and the evaporation process continues to occur at the same rate as in the pinned contact line case.Substituting (39) into (40), and using the dimensionless version of (4), which gives the volume loss due to the evaporation, we get the following equation for R n+1 , Employing the assumption that y n+1 = 1, we get the cubic equation which was prescribed in (12).

A7 Numerical sensitivity analysis
We have checked that the three parameters: roughness denoted by c, the slope of ln(l n ) versus n, where l n denotes the individual distances between rings, and the scaled position of the outermost deposit denoted by P, whose trends we present in Section 4, are not sensitive to the specific choices of the numerical parameters Dr, Dt, and Dr*, as far as convergence of the numerical scheme is achieved.To be more specific, we define convergence for the results presented in Section 4 under the requirement that the general characterisation of the deposited patterns changes by less than 10% for the pinned contact line case and the stick-slip motion of the contact line case and, for convenience, it changes by less than 20% for the unpinned contact line case following a decrease in Dr and Dt or Dr* by a factor of 2. The majority of the results, however, change by a few percent at the most upon the above modifications of the numerical parameters.For several numerical examples, see Tables 1-8.
Table 1 Sensitivity study in the case of stick-slip motion of the contact line with D = 0.25, y r = 0.4, and b = 10 (see the corresponding point on Fig. 9)

Fig. 1
Fig. 1 Cross-sectional illustration of an evaporating drop of a dilute solution on a solid substrate (not to scale), where we term the region between r = r* and the time-dependent position of the contact line r = R, of length Dr* R À r*, such that Dr*/R { 1, the vicinity of the contact line.

Fig. 2
Fig. 2 An enlarged cross-sectional view of the vicinity of the contact line which is positioned at R(t): (a) before the deposition and (b) after the deposition.The sketch visualizes our assumption that (hf)(r,t) = (hf)(r*,t) in the vicinity of the contact line, namely for r* r r r R(t), and illustrates the meaning of M(r,t) which is assumed to be flat in the region of length Dr*, where 0 o Dr*/R(t) { 1, and of the smooth function j(r,t) defined for 0 r r r R(t).

Fig. 3
Fig. 3 An example of the radial (r) variations of the deposited solute mass per unit area of the substrate M(r,T).

Fig. 4
Fig. 4 Examples for the deposited solute mass, M(r,T), subject to pinned contact lines for different relative roughness, c, levels.

Fig. 5
Fig. 5 Examples for the deposited solute mass, M(r,T), subject to unpinned contact lines for different relative roughness, c, levels.Note that the boundaries r n * of the domains for which we solve the advection-diffusion equation satisfy r min r r N * o r NÀ1 * o Á Á Á o r 1 * o r 0 *, where the simulation was stopped at r min = 0.15, so originally in the unpinned contact line case we have obtained a piecewise-defined function M(r,T).Thus, in order to obtain a smooth deposit, we have smoothed our solution by interpolating M(r,T) between r n * and r n+1 *, for n = 0, 1,. ..,NÀ 1.

Fig. 6 (
Fig. 6 (a and b) Examples for the deposited solute mass, M(r,T), subject to a stick-slip motion of the contact line for different levels of relative roughness, c, and of the position of the outer-most deposit, P, where in the inset we give the full results of the calculation and in (c and d) we give the respective ((c) for (a) and (d) for (b)) values of ln(l n ) n and the approximate slopes of these functions.

Fig. 7
Fig. 7 Roughness c versus the diffusion D for (a) a pinned contact line and (b) an unpinned contact line.The solid lines are to guide the eye along the calculated points.Note that in the unpinned contact line case, there is a missing point for b = 30 and D = 5; this is because the numerical scheme (for the unpinned contact line case) did not converge for these parametric values.

Fig. 8
Fig. 8 (a) Roughness c and (b) the position of the outermost deposit P versus the non-dimensional diffusion coefficient D for a stick-slip motion of the contact line, where y r = 0.6.The solid lines are to guide the eye along the calculated points.

Fig. 9
Fig. 9 (a) Roughness c versus the receding contact angle y r for a stick-slip motion of the contact line, where D = 0.25, y r = 0 corresponds to the pinned contact line case, and y r = 1 corresponds to the un-pinned contact line case, where we have smoothed the deposited shape M(r,T), and (b) (top) the slope of ln(l n ) as a function of n versus the receding contact angle y r and (bottom) the natural logarithm of the number of deposits ln(N) versus the receding contact angle y r .The solid lines are to guide the eye along the calculated points.

Table 3 Table 4
study in the case of stick-slip motion of the contact line with D = 0.25, y r = 0.95, and b = 30 (see the corresponding point on Fig. Sensitivity study in the case of stick-slip motion of the contact line with D = 5, y r = 0.6, and b = 20 (see the corresponding point on Fig. 8) Sensitivity study in the case of stick-slip motion of the contact line with D = 0.1, y r = 0.6, and b = 5 (see the corresponding point on Fig.
r o Rg r :

Table 5
8) Sensitivity study in the pinned contact line case with D = 0.1 and b = 15 (see the corresponding point on Fig. 7(a))

Table 6
Sensitivity study in the pinned contact line case with D = 5 and b = 25 (see the corresponding point on Fig. 7(a))

Table 7
Sensitivity study in the unpinned contact line case with D = 5 and b = 10 (see the corresponding point on Fig. 7(b))

Table 8
Sensitivity study in the unpinned contact line case with D = 0.1 and b = 5 (see the corresponding point on Fig. 7(b))