Sessile drop evaporation in a gap – crossover between diffusion-limited and phase transition-limited regime

Abstract We consider the time evolution of a sessile drop of volatile partially wetting liquid on a rigid solid substrate. The drop evaporates under strong confinement, namely, it sits on one of the two parallel plates that form a narrow gap. First, we develop an efficient mesoscopic long-wave description in gradient dynamics form. It couples the diffusive dynamics of the vertically averaged vapour density in the narrow gap to an evolution equation for the profile of the volatile drop. The underlying free energy functional incorporates wetting, interface and bulk energies of the liquid and gas entropy. The model allows us to investigate the transition between diffusion-limited and phase transition-limited evaporation for shallow droplets. Its gradient dynamics character allows for a long-wave as well as a full-curvature formulation. Second, we compare results obtained with the mesoscopic long-wave model to corresponding direct numerical simulations solving the Stokes equation for the drop coupled to the diffusion equation for the vapour as well as to selected experiments. In passing, we discuss the influence of contact line pinning.


Introduction
The dynamics of the liquid-vapour phase change, i.e., evaporation and condensation, plays a very important role in many systems involving films or drops of simple or complex liquids on solid substrates (Brutin 2015).Examples of practical importance include printing, coating and deposition processes (Brinker et al. 1992;Routh 2013;Thiele 2014), as well as cooling, moisture capturing and heat exchange technologies (Oron et al. 1997;Nguyen et al. 2018;Jarimi et al. 2020).In consequence, the evaporation of sessile drops of volatile liquids on rigid solid substrates is extensively studied in experiment and theory (Craster & Matar 2009;Cazabat & Guena 2010;Hu & Larson 2002;Semenov et al. 2011b;Erbil 2012;Kovalchuk et al. 2014;Larson 2014).
Thereby, the dynamics of droplet evaporation is controlled by the intricate interplay of various transport processes, namely, of heat and material within and between the liquid and the gas phase.They influence interface, temperature and concentration profiles, in turn causing pressure gradients as well as thermal and solutal Marangoni forces (Nepomnyashchy et al. 2002).These then drive convective motion within the liquid.For droplets on solid substrates, wettability and its interplay with evaporation in the region of the three-phase contact line also plays a crucial role (Plawsky et al. 2008).Although the involved processes can be modelled employing the full hydrodynamic description based on (Navier-)Stokes equations for the liquid and (advection-)diffusion equations for solutes in the liquid and vapour in the gas phase (Petsi & Burganos 2008;Bhardwaj et al. 2009), in many cases reduced descriptions are used.Common examples are long-wave (or lubrication, or thin-film) models for the liquid that are valid for small contact angles and interface slopes (Oron et al. 1997;Craster & Matar 2009;Saxton et al. 2016;Ji & Witelski 2018).
In all cases, the description of the dynamics of evaporating liquid films and drops on solid substrates crucially depends on the model for the evaporation rate j ev .It enters the kinematic boundary condition employed at the free liquid-vapour interface (Levich 1962;Leal 2007) and gives the mass loss per time and interface area.The rate j ev depends on material properties, thermodynamic state and on interface and system geometry (Oron et al. 1997;Plawsky et al. 2008;Erbil 2012).
One may distinguish two main approaches to the determination of j ev depending on the character of the process that limits the mass transfer across the interface.The limiting step can be either the actual phase change, e.g., for evaporation the transition of molecules from liquid state to gas state, or the diffusive transport of the vapour within the gas surrounding the drop, e.g., for evaporation transport away from the interface (Picknett & Bexon 1977;Sultan et al. 2004).Here, we call the two approaches (phase) transitionlimited and diffusion-limited, respectively.Other important distinctions are (i) whether the process is considered under homogeneously isothermal conditions or whether latent heat and heat transport are incorporated as further rate-limiting influences, and (ii) whether the evaporation is into pure vapour or into an inert gas.Only in the latter case one considers mass diffusion.
Diffusion-limited evaporation is considered in many models for evaporating liquid drops and films on solid substrates, either for simple liquids or suspensions and solutions.Such models are used and analysed, e.g., by Bourges-Monnier & Shanahan (1995); Deegan et al. (1997Deegan et al. ( , 2000)); Hu & Larson (2002); Cachile et al. (2002); Poulard et al. (2003); Erbil et al. (2002); Sultan et al. (2004); Popov (2005); Hu & Larson (2005); Shahidzadeh-Bonn et al. (2006); Murisic & Kondic (2008); Eggers & Pismen (2010); Semenov et al. (2011a); Tsoumpas et al. (2015); Saxton et al. (2016).An overview of earlier work is given by Hu & Larson (2002).This approach assumes that the phase transition is much faster than diffusion, i.e., directly at the liquid-gas interface the vapour is at saturation (Maxwell 1890;Langmuir 1918).In consequence, the local evaporation rate along the liquid-gas interface is controlled by the vapour diffusion within the entire gas phase.For shallow macroscopic drops, i.e., in the limit of small contact angles, the evaporation rate has a square-root divergence at the three-phase contact line (Deegan et al. 2000;Hu & Larson 2002;Popov 2005).† Note that, in contrast, drops with contact angle π/2 show a uniform rate, as can easily be seen from the electrostatic equivalent using "mirror charges".The effect that evaporative cooling has on the saturation concentration and the dependence of the diffusion constant on pressure is incorporated in Refs.(Dunn et al. 2008;Sefiane et al.Sessile drop evaporation in a gap 3 2009; Dunn et al. 2009b,a).Evaporating thin liquid films in a gap geometry are considered by Sultan et al. (2004Sultan et al. ( , 2005) ) -in the case of weak surface modulations a nonlocal single thin-film equation of closed form is determined employing Hilbert transforms.The influence of wettability and capillarity for mesoscale drops was considered by Eggers & Pismen (2010).There, a single, though non-local, equation for the dynamics of the thickness profile is obtained.A comparison of the diffusive and evaporative time scales was discussed by Ledesma-Aguilar et al. (2014) in terms of a Lattice-Boltzmann model of a volatile drop, yet assuming a diffusion slower than phase change.
The transition-limited case is considered in a number of different flavours: The "kinetic" approach by Burelbach et al. (1988); Joo et al. (1991) assumes a uniform constant saturated vapour density in the gas and determines the strength of evaporation/condensation via the difference of film surface temperature and the uniform saturation temperature in the gas phase.The approach is also normally applied if evaporation is into a pure vapour atmosphere, i.e., any vapour dynamics is then neglected.The derivation is based on a discussion of mass, energy and momentum flows across the liquid-vapour interface resulting, e.g., in the incorporation of vapour recoil effects.The approach is adopted in many later works, e.g., Anderson & Davis (1995); Hocking (1995); Oron & Bankoff (1999); Warner et al. (2003); Gotkis et al. (2006); Murisic & Kondic (2008); Savva et al. (2017).Dependencies of evaporation rate on interface curvature and wettability are normally not incorporated.
However, these effects are included in another flavour of the transition-limited case as presented by Potash & Wayner (1972); Moosman & Homsy (1980); Wayner (1993); Sharma (1998); Padmakar et al. (1999); Kargupta et al. (2001); Ajaev & Homsy (2006); Ji & Witelski (2018).It seems the earliest model for an evaporating meniscus influenced by Laplace (curvature) and Derjaguin (or disjoining) pressures is given by Potash & Wayner (1972), where evaporation into pure vapour is considered.At the liquid-vapour interface the vapour is at saturation and then varies vertically due to hydrostatic influences, i.e., it is always at equilibrium.The saturation pressure depends on Laplace and Derjaguin pressure.In consequence, Laplace pressure γκ (the product of liquid-gas interface tension γ and the interface curvature κ) and Derjaguin pressure Π enter the evaporation flux j ev , however, as argument of an exponential (Wayner 1993;Sharma 1998;Padmakar et al. 1999;Kargupta et al. 2001).Ultimately, evaporation is driven by a temperature difference between liquid and vapour.So, the process is seen as "transition-limited", but what limits the mass transfer between the phases is the diffusion of heat within the liquid.The gas phase itself is always at uniform constant temperature and pressure.Somewhat similar expressions are derived and/or used by Samid-Merzel et al. (1998); Ajaev & Homsy (2001);Lyushnin et al. (2002); Pismen (2004); Leizerson et al. (2004); Ajaev (2005a,b); Thiele (2010); Rednikov & Colinet (2010), to study dewetting volatile films, vapour bubbles in microchannels and evaporation fronts.There, however, the transitionlimitation is indeed due to the phase transition at the interface and the vapour is not at saturation, but at fixed vapour pressure (chemical potential).Furthermore, a direct proportionality of j ev and the sum of γκ and Π is used in the evaporation term.
The majority of works on thin-film models that incorporate evaporation either exclusively consider the transition-limited or the diffusion-limited case.A small number of works exist that compare the two approaches (Murisic & Kondic 2008, 2011).A partial comparison is also done for a model based on a Stokes description of the liquid (Petsi & Burganos 2008).The only work we are aware of that develops a more general model containing both limiting cases is by Sultan et al. (2005).
As laid out above, most evaporation models assume either (i) limitation by vapour diffusion in the gas phase or (ii) limitation by heat diffusion in the liquid phase, or (iii) The height-averaged density between the liquid and the upper plate is ρ(r, t) for the vapour, and ρair(r, t) for the remaining (inert) gas, i.e., the total height-averaged gas density is ρtot = ρ+ρair.
There can be exchange between the liquid and the vapour due to evaporation/condensation.
limitation by mass transfer between liquid and gas phase.Case (i) is not applicable for evaporation into pure vapour and implicitly assumes uniform total pressure in the gas phase.Any pressure gradient would trigger convective flows in the gas phase what is, however, excluded.Wetting-and capillarity-influence on saturation concentration can be incorporated.Case (ii) assumes a uniform vapour concentration corresponding to the value at saturation at a gas reference temperature that differs from the temperature of the liquid at the liquid-gas interface.This difference drives evaporation.Wettingand capillarity-influence on evaporation can be incorporated.Finally, case (iii) assumes constant vapour pressure (chemical potential) in the gas phase what results in inhomogeneous evaporation due to wetting-and capillarity-dependencies.Any evaporation-induced inhomogeneous vapour and total gas pressure are assumed to instantaneously equilibrate.Nearly all mentioned models focus on one of these cases and do not allow for an analysis of transitions between the different limiting cases.Our present aim is to develop a relatively simple long-wave model in gradient dynamics form that bridges cases (i) and (iii) for the specific geometry of a sessile drop of partially wetting liquid evaporating into the narrow gap between two parallel rigid smooth solid plates (see Fig. 1).There, the coupled liquid and gas dynamics can be described by kinetic equations of reduced dimensions.For simplicity, we consider a completely isothermal system -thermal effects can be incorporated later on.The simple gradient dynamics approach shall allow one to incorporate our model as a building block into a wide class of thin-film models for more complex settings.Furthermore, we show that the developed simplified mesoscopic approach favourably compares to a full macroscopic description as well as to experiments.
As a macroscopic model we employ a Stokes model coupled to diffusion in the gas phase.Thereby, an evaporation rate j evap is implemented into the boundary condition at the liquid-gas interface that is equivalent to the one used in our gradient dynamics model.At the contact line, a Navier slip condition is used.In the gas phase a vapour diffusion model is employed that fully resolves the space within the gap.Versatile variants of this model have been successfully used to account for e.g.multi-component droplet evaporation (Diddens et al. 2017;Li et al. 2019) or droplets evaporating on a thin oil film (Li et al. 2020).
Experimental results on evaporating sessile droplets are extensive (see the reviews of by Cazabat & Guena (2010); Brutin & Starov (2018); Zang et al. (2019) ), but only a limited number of the previous studies investigated the effect of confinement inside microfluidic channels (Bansal et al. 2017a,b;Hatte et al. 2019) or in a box (Shahidzadeh-Bonn et al. 2006).To provide counterpart experiments for our theory, we analyse the evaporation of a droplet between two horizontal plates, where confinement is only imposed in one (vertical) direction.
This paper is structured as follows: In Section 2 we present the theoretical and experimental approaches that are compared in this work.In particular, Section 2.1 discusses the general form of long-wave gradient dynamics models for one and two scalar fields with combined conserved and non-conserved dynamics.Then, in Section 2.2 we derive the gradient dynamics model for the evaporating drop in the considered small-gap geometry, and specify all necessary parameters and specific functions.The subsequent Sections 2.3 and 2.4 briefly introduce Stokes description and experimental setup, respectively.Next, Section 3 presents results obtained with the developed thinfilm model that are in the subsequent Section 4 compared with Stokes-equation results and corresponding experiments.Finally, Section 5 concludes with a discussion of the limitations of the presented approach and an outlook toward its further development and application.

Long-wave gradient dynamics models
The dynamics of a layer or shallow drop of nonvolatile liquid in long-wave approximation (Oron et al. 1997;Craster & Matar 2009) is characterised by the evolution of a single field -the layer thickness.As first noted by Oron & Rosenau (1992) and Mitlin (1993), the corresponding dynamic equation for the layer thickness h can be written in gradient dynamics form for a conserved field (see appendix of Thiele (2018) for a derivation of this form via Onsager's variational principle (Doi 2011)).To account for evaporation in the phase transition-limited case (Lyushnin et al. 2002), one adds thermodynamically consistent nonconserved contributions to the dynamics (Thiele 2010) and obtains where the energy functional F[h] contains wetting energy and surface energy of the free liquid-gas interface.In the simplest case, p vap is the imposed constant external vapour pressure in the gas phase (also the total gas pressure is then constant if evaporation is not into pure vapour).Here and in the following ∂ t denotes the partial time derivative and ∇ = (∂ x , ∂ y ) T is the two-dimensional (2d) spatial gradient operator.The functions Q(h) 0 and M (h) 0 are the positive mobilities of the conserved and the non-conserved contributions, respectively.For a discussion of different forms of M (h) see Thiele (2014).
For systems with more degrees of freedom, the described one-field model (2.1) is extended by incorporating the dynamics of further fields (Thiele 2018).In the context of thin-film hydrodynamics, two-field gradient dynamics models are presented and analysed for (i) dewetting two-layer films on solid substrates, i.e., staggered layers of two immiscible fluids (Pototsky et al. 2004(Pototsky et al. , 2005;;Jachalski et al. 2013;Bommer et al. 2013), (ii) decomposing and dewetting films of a binary liquid mixture (with non-surface active components) (Thiele 2011;Thiele et al. 2013;Diez et al. 2021), (iii) the dynamics of a liquid film that is covered by an insoluble surfactant (Thiele et al. 2012(Thiele et al. , 2016)), (iv) the spreading of a liquid drop on a polymer brush (Thiele & Hartmann 2020) and (v) on an S. Hartmann, C. Diddens, M. Jalaal and U. Thiele elastic substrate without (Henkel et al. 2021) and with (Henkel et al. 2022) Shuttleworth effect.In all these cases, the model is of the form where the indices a, b = 1, 2 refer to the two fields.For the considered relaxational dynamics, the Q(u 1 , u 2 ) and M(u 1 , u 2 ) represent 2 × 2 positive definite and symmetric mobility matrices for the conserved and non-conserved parts of the dynamics, respectively, written here in terms of their components Q ab and M ab .In examples (i) to (iii), both fields show a conserved dynamics, i.e., M = 0.However, in general, they can also show a purely nonconserved dynamics, i.e., Q = 0 (corresponding to a two-field model-A in the Hohenberg & Halperin (1977) classification), or a mixed dynamics as in cases (iv) and (v), i.e., M, Q = 0.The mobilities Q enter the fluxes j a = − 2 b=1 Q ab ∇(δF/δu b ) of the conserved part of the dynamics for both fields u a .They are given as linear combinations of the influences of both thermodynamic forces −∇(δF/δu b ).The components of M give the transition rates between the two fields and between the fields and the surroundings.The conserved fields u 1 and u 2 represent in case (i) the lower layer thickness h 1 and overall thickness h 2 , respectively (Pototsky et al. 2004(Pototsky et al. , 2005;;Jachalski et al. 2013) or the lower and upper layer thickness (Bommer et al. 2013).In case (ii), u 1 and u 2 represent the film height h and the effective solute height ψ = ch, respectively, where c is the height-averaged solute concentration, while in case (iii), u 1 and u 2 represent the film height h and the surfactant coverage.respectively(Thiele et al. 2012).Finally, in cases (iv) and (v) u 1 represents the drop height while u 2 stands for the local amount of liquid in the polymer brush (Thiele & Hartmann 2020) and the elastic-liquid interface profile (Henkel et al. 2021), respectively.2.2.Gradient dynamics for volatile liquid in small-gap geometry

Gradient dynamics form
Having set the stage for formulating thin-film models as gradient dynamics on an underlying energy functional, we next introduce such a thin-film model for an evaporating sessile liquid drop with profile h(r, t) in a gap of width d (see Fig. 1).We employ two fields, namely, on the one hand, the amount ψ 1 (r, t) of the substance in liquid state in the drop per substrate area and, on the other hand, the amount ψ 2 (r, t) of the substance in vapour form in the gas phase also per substrate area.The field ψ 1 (r, t) is proportional to the thickness of the liquid film: where ρ liq is the constant liquid density (to be specified later).All employed densities are number densities, i.e., are given in units of particles per volume.The field ψ 2 (r, t) is proportional to the height-averaged vapour density ρ(r, t), namely, (2.4) The gas phase can either consist of pure vapour or of vapour in an inert gas (called "air" in the following).The latter has height-averaged density ρ air (r, t) and the resulting total gas density is ρ tot = ρ + ρ air .The literature distinguishes three main cases: (i) Evaporation into pure vapour treated isothermally.In this case no diffusion occurs, all dynamics in the gas phase is due to pressure equilibration via convective motion.However, this is a very fast process as it occurs with the speed of sound (Maxwell 1890).
On the time scale of evaporation one then assumes uniform vapour, i.e., gas pressure.This implies ρ(r, t) is uniform.
(ii) Evaporation into air treated isothermally.In this case vapour diffusion is very important.As in case (i), the total pressure equilibrates fast, i.e., also ρ tot is constant and uniform, e.g., p = k B T ρ tot for an ideal gas.In consequence, (2.5) (iii) As case (i) or case (ii), but not treated isothermally.Then, heat diffusion becomes a limiting factor as the heat used up as latent heat during the liquid-vapour phase transition needs to be transported to the interface.It also becomes important, that normally a jump in the heat flux is considered at the interface that ultimately controls the evaporation flux.Here, we will not discuss this case.Next, we write the model for the coupled dynamics of the local amounts of liquid and vapour in the gradient dynamics form (2.2).In particular, the ψ i follow the mixed conserved and nonconserved dynamics Note, however, that h is the relevant field for many of the physical effects and that often it may be more convenient to write the governing equations in terms of h [Eq.( 2. 3)] and ψ = ψ 2 .First, we focus on case (ii) before discussing the amendments needed for case (i).As on the considered time scales, the total gas pressure is uniform, there is no pressure gradient that would drive convective transport in the gas phase.In consequence, there is no dynamic coupling between liquid and gas layer, i.e., Q 12 = Q 21 = 0. Transport in the gas phase is then limited to vapour diffusion within the air.As a result, the conserved mobilities are † where η is the dynamic viscosity of the liquid and D is a diffusive mobility constant.The latter is related to the usual diffusion constant D of the vapour particles in air by where k B is Boltzmann's constant and T the temperature.In the simplest case, the nonconserved mobilities are where M is an evaporation rate constant, which can be estimated e.g. from the Hertz-Knudsen equation (Knudsen 1915;Librovich et al. 2017).The particular matrix ensures that phase change is driven by differences in chemical potentials and that ψ 1 + ψ 2 is conserved, i.e., their sum satisfies a continuity equation The form of the diffusive mobility is for solutions discussed by (Thiele 2011;Xu et al. 2015).The form Dψ2 = D(d − h)ρ is the direct equivalent in the present case.
S. Hartmann, C. Diddens, M. Jalaal and U. Thiele The free energy in long-wave approximation is where the f i 's are bulk liquid, vapour and air energies per volume, γ is the liquid-gas interface tension and g(h) is the wetting energy per area.The functional is accompanied by the constraint of particle number conservation across the two phases.The condition is where n is the mean (liquid and vapour) particle number per substrate area.Here, particle flux through the boundaries of the domain Ω is excluded, however, it can be easily incorporated.Also, gravity is neglected as we consider small droplets, but may be added in the form of potential energy.
The presented long-wave formulation of the mesoscopic model is best suited for shallow droplets.However, using a full-curvature trick similar to (Gauglitz & Radke 1988;Snoeijer et al. 2007) one may also study the evaporation of droplets with larger contact angles.Then, the surface energy term in Eq. (2.9) is written with the full metric factor γ 1 + (∇h) 2 and the evaporation rate per surface area M in Eq. (2.8) is replaced by M 1 + (∇h) 2 .Below we will refer to this amendment as the full-curvature formulation of the mesoscopic model.
Further we remark that Eq. (2.9) is in a somewhat mixed form as it considers at the same time equations of state encoded in the f 's for liquid and vapour that may result in phase change, but also uses the film height h as it is the most important quantity for mesoscopic hydrodynamics.As stated above, we treat ρ liq as constant and will consider the vapour as ideal gas.This is the result of a two-step procedure, explained in the next section (that may be skipped by the reader who is mainly interested in the resulting thin-film equations).

From real to ideal vapour
An equation of state for a real gas, e.g., a van der Waals gas, predicts coexisting liquid and gas densities at equilibrium.In an out of equilibrium hydrodynamic description, these densities may then vary in space and time.However, as there is no simple way to express the three fields ρ liq (r, t), ρ(r, t) and h(r, t) in terms of the two fields ψ 1 (r, t) and ψ 2 (r, t) we follow a two-step procedure: (i) We assume a thick homogeneous sharp-interface flat film of thickness h and calculate ρ liq and ρ vap at coexistence (as functions of temperature T ) using Using the trick to treat ξ = d − h as independent field, we minimise with respect to ρ liq , ρ vap , h and ξ.Here, μ and p are Lagrange multipliers for mass and volume conservation.We obtain (2.13) i.e., the standard Maxwell construction for phase coexistence.
In this step we can either employ a single function f = f vap = f liq that allows for a liquid-gas phase transition (e.g. the van der Waals free energy) or use a purely entropic f vap (ρ) = k B T ρ[log(Λ 3 ρ) − 1] and combine it with a f liq that allows for a phase transition.These we use to obtain the coexisting ρ liq and ρ vap analytically or (most likely) numerically.From hereon, symbols ρ liq , ρ vap and f liq denote the values at coexistence.
(ii) Next we approximate the equation of state in the vapour and the liquid phase.For the liquid phase we neglect any compressibility, i.e., we fix the density in the liquid layer to ρ liq .With other words the "liquid branch" of the equation of state p(ρ) is replaced by a vertical line at ρ = ρ liq .The "gas branch" of the equation of state is either directly replaced by the ideal gas law given in the previous paragraph or, alternatively, is expanded up to linear order about ρ = ρ vap .This results in a shifted and scaled ideal gas law.The latter approach has the advantage that coexistence pressure and concentration are exactly as in the equation of state one started with.
In this way the relation h = ψ 1 /ρ liq introduced above becomes meaningful and the relation of variations w.r.t.h and ψ 1 alluded to earlier is justified.Also using ρ = ψ 2 /(d− h), the free energy functional Eq.(2.9) is written as Here, it only depends on ψ 1 and ψ 2 .Note that f liq is now a constant given by the Maxwell construction.Next we bring all the information together and present the thinfilm equations.

Thin-film equations
The energy functional (2.14) in terms of ψ 1 and ψ 2 is now minimised together with the particle number constraint (2.10) (Lagrange multiplier µ) w.r.t.variations in the two fields.This gives where the brackets contain Laplace pressure, Derjaguin pressure, liquid energy and vapour pressure.Note that the somewhat unusual final contribution to the vapour pressure is a direct consequence of the assumption of constant p tot .Here and in the following, we use h and ρ as abbreviations where appropriate.The variation w.r.t.ψ 2 gives i.e., a difference in chemical potentials.If we now specify vapour and air to be ideal gases, with the mean free path S. Hartmann, C. Diddens, M. Jalaal and U. Thiele length Λ, we obtain The gradient of this pressure drives the dynamics in the liquid film.The constant liquid energy density f liq and the constant ideal pressure k B T ρ tot do not contribute to the conserved dynamics, but the final term in the curly parenthesis does.It is a direct consequence of treating case (ii), i.e., of imposing a uniform gas pressure, i.e., constant ρ tot .The other variation is then (2.18) Also here the second term in the curly parenthesis is a consequence of the imposed uniform pressure in case (ii).
Introducing the obtained expressions into the general two-field gradient dynamics (2.6) with (2.7) and (2.8) gives This is the final result for case (ii).
Then, the saturation vapour density ρ sat above a flat thick film is obtained by setting the transfer term E = 0 and dropping capillarity and wettability influences: The corresponding film height H is determined by the conservation of mass n = The mentioned somewhat unexpected terms in (2.19) are well behaved: The additional contributions in the equation for ψ 2 provide a factor 1/[1 − ρ/ρ tot ] to the generalised diffusion constant.Normally, ρ/ρ tot 1 and can be neglected.If, however, ρ/ρ tot → 1 we approach the limit of pure vapour where diffusion is not the proper transport process any more: the factor in question diverges, i.e., ρ becomes instantaneously uniform.The additional term in the equation for ψ 1 corresponds to a flux proportional to 1 1−ρ/ρtot ∇ρ.For ρ/ρ tot 1 this is the gradient in partial vapour pressure that drives some flow in the adjacent liquid layer, one can see it as an "osmotic coupling".It is normally very small as compared to the other pressure gradients.The local evaporation rate M E contains additional terms proportional to the ratio of total gas density and liquid density ρ tot /ρ liq that we expect to be small.For dry air and water the ratio is about 10 −3 , and humid air has an even lower density than dry air.
For ρ/ρ tot 1 and ρ tot /ρ liq 1, Eqs. (2.19) written in terms of h and ρ reduce to Here, we have introduced the diffusion constant D = Dk B T and the evaporation rate constant of M = M /ρ 2 liq , thus converting the local particle evaporation rate M E into a volume rate j evap = M E/ρ liq .Equations (2.21) allow one to study the crossover between transition-limited drop evaporation dynamics (small ρ) and diffusion-limited evaporation dynamics (ρ at drop close to saturation ρ sat ).
In case (ii) considered up to here, convective motion in the gas layer is neglected assuming a uniform total gas pressure p tot , i.e., a uniform total gas density ρ tot = ρ air + ρ that is an externally controlled parameter of the system.However, this directly couples ρ air to the vapour density ρ what itself results in small additional contributions to the pressure in the liquid, to the evaporation/condensation rate and vapour diffusion.All these terms are direct consequences of the gradient dynamics structure.Although one may neglect them due to their smallness one needs to keep in mind that this breaks the thermodynamic consistency and, therefore, can result in unphysical behaviour.
If we consider case (i), evaporation into pure vapour, diffusion is excluded right from the beginning.As in case (ii), we assume that convective motion is very much faster than evaporation, what directly implies that the vapour density ρ is uniform in the entire gas layer.It is set by the conditions at the lateral boundaries of the gap.With other words, the vapour pressure is then a given constant and the governing thin-film equation in gradient dynamics form is Eq.(2.1).The corresponding energy is Eq.(2.9) without the air energy and with constant ρ = ρ vap , i.e., Minimisation w.r.t.film thickness gives i.e., the thin-film equation (2.1) becomes Grouping all constants in the evaporation term into a single one, the equation corresponds to the standard thin-film description in the transition-limited case (Pismen 2004; Thiele 2010).

Specific functions and parameters
One may now proceed by nondimensionalising, thereby using sensible values of ρ liq and f liq as parameters.However, as they are not independent, this may result in artefacts.The better approach is to employ a specific equation of state / free energy that allows for liquid-gas phase transition and calculate these values for a particular temperature T .
An example is the van der Waals equation of state, see e.g., Landau & Lifshitz (1980, § 76 & 84) giving the pressure where a and b are constants (attraction strength and effective excluded volume, respectively).Approximated, in the gas phase with a, b → 0 it becomes the ideal gas law p id = k B T ρ.The corresponding Helmholtz free energy per volume is Further, the saturation vapour density ρ sat follows from the equilibrium condition, e.g., in the simple case of a thick flat liquid film (no Laplace and Derjaguin pressure) setting the evaporation rate to zero in Eq. (2.21): (2.27) The saturation vapour pressure p sat resulting from the equation of state may then be used to express the humidity as φ = p/p sat .
If one assumes a specific wetting potential g(h) that allows for partial wetting, a second spatially homogeneous equilibrium state is found typically at very small height h a that corresponds to a microscopic adsorption layer (2.28) The height h a depends on the vapour density ρ and will be slightly shifted from the height h p where the Derjaguin pressure vanishes (g (h p ) = 0), i.e., the adsorption layer height of the saturated case h p = h a (ρ sat ).In fact, this allows for a set of spatially homogeneous equilibrium states where the vapour density is an arbitrary constant 0 < ρ ρ sat controlled externally, e.g., by the humidity of the laboratory, while the substrate is macroscopically dry (the "moist case" of De Gennes (1985)) with h = h a .Note that due to the saturation condition Eqs.(2.27) and (2.28) we find g (h a ) 0, i.e., the adsorption layer height is shifted to lower values h a h p .When, on the other hand, an initial ρ(r, t) ρ sat is given, no equilibrium is guaranteed to exist (depending on the choice of the wetting potential).The evaporation rate E will then become negative leading to a transfer of material from the vapour phase to the liquid film, i.e., condensation.
In the following, we use the simple wetting potential for partially wetting liquids with the Hamaker constant H A ≈ 5 3 γh 2 p θ 2 e that relates to the macroscopic equilibrium contact angle θ e by a mesoscopic Young's law: cos θ e = 1 + 1 γ g(h p ) (Churaev 1995).For the simulation of the thin-film model in section 3 we use the following set of realistic parameters (unless stated otherwise): initial drop height h 0 = h(r = 0, t = 0) = 1 mm, gap height d = 3 mm, precursor film thickness h p = h 0 /1000, diffusion constant D = Dk B T = 2.8 × 10 −5 m 2 /s (Lee & Wilke 1954;Marrero & Mason 1972), temperature T = 22 °C, contact angle θ e = 44.38°,and a domain radius of L = 25 mm.The liquid properties are chosen according to water at room temperature: viscosity η = 8.9 × 10 −4 Pa s, surface tension γ = 72 × 10 −3 N/m, particle density according to mass Sessile drop evaporation in a gap 13 density and molar mass ρ liq,mass = 997 kg/m 3 18 kg/mol N A , and the saturation vapour pressure p sat = 2643 Pa.The parameter values are extracted from Lide (2004).The constant of phase change M is employed as a free parameter that allows us to move between the different limiting cases.

Radial symmetry and numerical implementation
For an efficient numerical simulation of the dynamics, we consider a radial symmetry, i.e., for the long-wave model h(x, t) = h(r, t) and ψ(x, t) = ψ(r, t).This allows us to reduce equations (2.21) to a spatially one-dimensional model in polar coordinates where The dynamic equations are then solved using the finite-element method implemented in the C++ library oomph-lib (Heil & Hazel 2006), which employs a second order BDF scheme for temporal integration.In particular, the library provides both spatial and temporal adaptivity, which is essential due to the multi-scale character of the dynamics, e.g., when spatially resolving the fields in the bulk and near the contact line.In the simulations, the domain is typically discretised with ≈ 1000 grid points with spacings that range from 2 −7 to 2 −18 times the domain size.
In these calculations, we employ an adsorption layer of height h p that is by a factor of 10 3 smaller than the initial drop height h(r = 0, t = 0).Larger ratios are possible but result in a strongly increasing numerical effort.As the effective adsorption layer height is coupled to the gas phase [see e.g.Eq. (2.28)], the overall pressure balance can cause liquid transport through the layer during equilibration.The resulting fluxes are negligible if the adsorption layer is very thin compared to the droplet size.The effect can be further suppressed by modulating the evaporation rate coefficient M (h) with a (smooth) step function, effectively disabling evaporation from the adsorption layer.This is in particular important on the slow timescale of droplet evaporation.
To ensure smoothness of all fields at the drop center (r = 0), we employ natural (homogeneous) Neumann boundary conditions (BC) (2.32) In particular, these conditions allow for a nonzero vapour flux through the outer boundary into the laboratory environment.The external vapour concentration is here chosen such that it corresponds to a low relative humidity of ρ lab /ρ sat = 10 %.
S. Hartmann, C. Diddens, M. Jalaal and U. Thiele Results obtained with the developed thin-film model are presented in Section 3. In the subsequent Section 4 they are compared with Stokes-equation results and corresponding experiments.

Stokes description
Next, we develop a description in terms of the Stokes equation for volatile liquids in the gap geometry again taking into account diffusion processes in the gas phase and the mass-transfer at the liquid-gas interface.This will allow us to study the transition between transfer-limited and diffusion-limited cases also in the macroscale model.Then, important cases are compared to the thin-film model developed above.
The Stokes equations are solved on an axisymmetric domain, i.e., for the velocity field u(r, z, t) and the pressure p(r, z, t), we solve (2.33) in the liquid domain.At r = 0 the conventional boundary conditions u r = 0 and ∂ r p = 0 are imposed.There is no precursor film considered in the Stokes description.Instead, a Navier slip boundary condition is used at the substrate at z = 0.For comparison, the slip length is chosen to coincide with the precursor film thickness of the corresponding thin-film simulations, i.e., see Savva & Kalliadasis (2011) for a comparison of precursor film and slip length models.The liquid-gas interface is not described by a height function h(r, t), but by a parametric representation R(s, t).On this interface, the surface tension is imposed as normal traction, i.e., n • −p1 + η ∇u + (∇u) t = (−γκ) n , (2.36) where κ = −∇ • n is the full curvature and n is the normal vector.The kinematic boundary condition has to consider the evaporation rate j evap , i.e. the relative normal motion between fluid velocity u and interface motion Ṙ is enforced to be (2.37) The evaporation rate reads here where the local pressure p within the liquid appears instead of the terms −γ∇ 2 h + g (h).
In contrast to the thin-film model, where the Derjaguin pressure contribution g (h) can strongly influence the evaporation rate within the transition region between drop and precursor film, here, the singularity of the evaporation rate near the contact line is only limited by the finite mobility M in the Stokes description with slip-length.
The contact angle θ is weakly imposed at the contact line by employing γ(cos θ e r − sin θ e z ) as surface tension force at the contact line, which balances when the free surface attains the prescribed contact angle θ.
The vapour diffusion in the gas phase is also fully resolved in the direction orthogonal to the substrate.Thus, instead of solving for the evolution of the height-averaged vapour density ρ(r, t), the diffusion equation is solved for the density ρ(r, z, t).
No-flux conditions are used at the top plate, at the axis of symmetry and at the substrate-gas interface, which ranges from the contact line at (r c , 0) to the end of the domain at (L, 0), and lab humidity is imposed at the far end, i.e. (2.40) Finally, at the liquid-gas interface, the evaporation flux is considered: (2.41) The equations are implemented in oomph-lib employing a sharp-interface Arbitrary Lagrangian-Eulerian method, i.e. the mesh nodes are moving in such a way that the liquid-gas interface is always conforming with the mesh.Mesh reconstruction and subsequent interpolation of all fields to the new meshes is invoked whenever the mesh distortion exceeds given thresholds.The mesh has enhanced resolution at the contact line and the typical number of degrees of freedom of the discretised system is chosen to be around 1,000,000.

Experiments
Experiments are done using a simple shadowgraphy setup.The millimetric water droplet (volume of ≈ 0.4 µl) is placed on a cover glass with a thickness of ≈ 0.2 mm (initial contact angle of ≈ 65°) -see Fig. 2. The cover glass is placed on top of a holder such that the air occupies the space below it.The top surface is a thick plastic circular disk with a diameter (w) much larger than the size of the droplet.The gap between the plates (d) is accurately controlled with a micrometer positioning stage, and special care is taken to make sure both top and bottom surfaces are horizontal.An LED light source illuminates the test section, and the images of the droplet are recorded using a CMOS camera attached to a zoom lens (see Fig. 2 for example snapshots).Experiments are performed for different levels of confinement, varying the gap height d.The changing shape of the droplet over time is acquired by the standard processing of the images.

Results of thin-film model
First, we employ the thin-film model in radially symmetric form (2.30) with BC (2.31) and (2.32) to simulate an evaporating droplet for different modi of evaporation.are fixed.All simulations were initialised with a sessile droplet profile of 1 mm height, representing a droplet equilibrated without evaporation, in a homogeneous atmosphere of a constant (low) humidity.The snapshots are taken after initial transients have passed and a (quasi-static) vapour concentration profile in the gas phase has been established that subsequently changes on a much slower time scale.
In particular, the bottom panels show the drop profile (solid black line, right hand side scale), the vapour concentration (i.e., relative humidity) as bluish background shading and also as black dotted line (left hand side scale).The upper panels give the corresponding evaporation rate profiles on a logarithmic scale.The left column [Fig.3 (a)] presents the diffusion-limited case at relatively large M = 4 × 10 −10 m/Pa s, i.e., the diffusive transport of vapour through the gas phase is much slower than the phase change, thereby effectively controlling the entire process.Note that above the droplet the vapour is nearly saturated due to the gap-geometry.Starting in the contact line region the concentration logarithmically decays towards the edge of the plates where the humidity is kept at 10%.The analytical form is discussed below in appendix A.
Remarkable is the clearly visible divergence of the evaporation rate when approaching the contact line from inside the drop [top panel of Fig. 3 (a)].In this region the local evaporation rate changes exponentially by about seven orders of magnitude leading to a very sharp spike.See Fig. 11 in the appendix for a highly zoomed version of the peak demonstrating the smoothness of the solution in the vicinity of the contact line.
In contrast, Fig. 3 (c), i.e., the right column, shows the transition-limited case at relatively small M = 4 × 10 −18 m/Pa s, i.e., the diffusive transport in the gas phase is much faster than the phase change.Now it is the latter that controls the entire process.As diffusion is fast, all surplus vapour is rapidly transported out of the gap between the plates and the vapour concentration profile remains nearly homogeneous at lab humidity.In consequence, the evaporation rate is nearly constant along the entire surface of the drop and rapidly falls to zero in the contact line region [top panel of Fig. 3 (c)].
Beside the two limiting cases respectively considered by the two groups of thin-film models in the literature, our model is also able to simulate 'mixed cases' where the time scales of the involved processes are not strongly separated.An example is presented in the centre column of Fig. 3 at moderate M = 4 × 10 −15 m/Pa s.Here, diffusion is fast enough to keep the concentration at the drop centre below saturation.It is also sufficiently slow for a decaying concentration profile to develop outside the drop.Accordingly, also the profile of the evaporation rate [top panel of Fig. 3 (b)] shows features of both limiting cases: the flux is moderately large above the drop, increases by about half an order of magnitude towards the contact line region where it steeply decays to zero.
After having established the stark difference in evaporation flux and concentration profiles seen for single snapshots in the two limiting cases, next we consider how the time evolution differs between them.To do so, in Fig. 4 we study how drop volume and height change in time for a number of different gap widths d in the diffusion-limited and the phase transition-limited case.In the former case [Fig.4 (a)], the drop volume shows a linear decrease in time.Thereby the rate is proportional to the gap width d as it directly controls the overall vapour flux in our height-averaged setting (cf.Appendix A).Correspondingly, the drop height shrinks with a power law ∼ (t evap − t) 1/3 as expected for a nearly constant contact angle.Here, t evap is the time when the complete drop has evaporated.
In contrast, in the transition-limited case [Fig.4 (c)], the drop height decreases linearly S. Hartmann, C. Diddens, M. Jalaal and U. Thiele with a rate nearly independent of the gap width.Correspondingly, the drop volume shrinks with a power law (t evap − t) 3 .Note that in the final stage when the drop is rather small, the remaining volume is rapidly absorbed into the adsorption layer.This is a mesoscale effect resulting from the dominance of the wetting energy for very small drop size.We do not further discuss it here.Interestingly, in the intermediate case of moderate evaporation [Fig. 4 (b)], features of both limiting cases are visible: in contrast to the transition-limited case, the evolution depends on gap width, but less so than seen in the diffusion-limited case.Neither the drop volume nor the drop height decrease linearly.Instead, in the course of the time evolution they cross over from a diffusion-limited start, where slopes of approximately linear decays in volume are proportional to gap width, to a transition-limited end with cubic decrease of volume and final drop collapse.
In the transition-limited case, the phase transfer rate decreases with the drop size, as it is limited by its surface area.Note that this is not true for nanoscopic droplets, where the Laplace pressure significantly contributes to the drop chemical potential, causing an increased evaporation, i.e., the Kelvin effect.This is not observed in our mesoscale simulations, due to the size of the droplets.Nevertheless, we expect that the behaviour of very small droplets should always tend towards the transition-limited behaviour because their surface area shrinks while the rate at which lateral diffusion transports vapour particles away remains (nearly) constant.The magnification in the insets of Fig. 4 (b) therefore depicts the measured volume and height data shifted in time such that the curves coincide at the time t evap of complete evaporation.There, in particular at large gap heights d, i.e., when the diffusion is less dominant, for small drop sizes the behaviour converges to a common curve that resembles the transition-limited case of Fig. 4 (c).
Note that in all considered cases evaporation is sufficiently slow to only see a minor evaporation-induced difference (Morris 2001;Todorova et al. 2012;Rednikov & Colinet 2017) below 2% between the observed quasi-steady contact angle and the equilibrium angle.Further, the contact line smoothly recedes as no contact line pinning occurs for the assumed ideally homogeneous and smooth substrate, i.e., there is no contact angle hysteresis.This aspect will be refined in the following two subsections.

Dependence on the contact angle
Before we investigate the effects of contact angle hysteresis, we check the overall influence of the contact angle on the dynamics.We therefore perform long-wave simulations of droplets with varying equilibrium contact angle, i.e., we control the parameter θ e in the wetting potential (Eq.(2.29)) and initialise the simulation with a droplet of according equilibrium shape.Droplets of identical volume but different radii should then show different evaporation rates merely due to their differences in surface area, which masks the actual effect of the contact angle.Hence, in the following we consider droplets of identical initial radius and accept that they have various initial volumes, depending on their contact angle.Naturally, a drop of smaller volume evaporates faster than a large drop, again hiding the mere influence of the contact angle.Yet, we can account for the different drop volumes by normalizing the time scale with the initial volume t = t/V (t = 0).
We then employ measures identical to the ones in Fig. 4, but vary the contact angle θ e instead of the gap height, which now remains fixed at d = 1.5 mm.All other parameters are kept at the same values as given at the end of Section 2.2.4.In Fig. 5 we plot the measured drop volume V (t) and height H(t) versus the normalised time t = t/V (t = 0) for three different values of the phase transition coefficient M , again accounting for the  Notably, in all configurations the dynamic behaviour is mostly independent of the contact angle.The only visible differences occur when the droplets become very small.This is most pronounced for low contact angles and in the final drop collapse stage of the transition-limited case.This is again a mesoscale effect due to the dominant wetting potential for droplets of low film heights.The otherwise coinciding curves in Fig. 5 demonstrate that all effects of the contact angle on the evaporation rate are due to changes in the liquid surface area.In particular, low contact angle droplets have a higher surface area than a droplet with higher contact angle of the same volume and will evaporate faster, depending on the relation of the total evaporation rate to the surface area.See the appendix A for an asymptotic approach to determining the latter.There, we find that the total evaporation rate scales as 1/R 2 with the drop radius in the phase transition-limited case and is mostly independent of it in the diffusion-limited case, i.e., the behaviour is also independent of the contact angle.

Effects of contact line pinning
In all simulations discussed up to here, the contact line continuously recedes during evaporation in order to retain the contact angle governed by the wetting potential g(h) [Eq. (2.29)].However, our choice of a simple wetting potential does not account for contact angle hysteresis, e.g., due to roughness of the solid substrate.Hence, the previous discussion does not include any effects of pinned contact lines.
Inhomogeneities of the wettability are, however, easily incorporated into the formalism  but with an initial pinning of the drop's contact line due to a contact angle hysteresis of 30°.For comparison, the simulations without pinning are included as thin dashed lines.The noticeable corner in the height curves indicate the depinning event.While the contact line is pinned, the volume evaporation rate dV /dt remains nearly constant, accelerating evaporation particularly in the transition-limited case.
through a spatial modulation of the wetting potential (Thiele et al. 2003).In the following, we simulate a contact angle hysteresis by modulating g(h) with a (possibly smoothed) step function ξ(r) that drastically increases the wettability (by lowering the energy at the minimum of the wetting potential) in the substrate region of the initial droplet radius R 0 .Namely, This corresponds to a contact angle hysteresis with the threshold angles θ min = √ 0.1 θ max , i.e., between 14.03°and 44.38°in our simulations.
In Fig. 6 we provide a comparison between the long-wave simulations without pinning of the contact line, namely, the results of Fig. 4, with simulations including contact angle hysteresis.Apparently, independently of the considered parameter regime, the drop volume initially decreases linearly, i.e., the evaporation rate appears constant.Note that here as well the drop height seems to decrease nearly linearly, because the volume V and height H of a shallow drop of fixed base radius R are related as V = π/2 HR 2 .At a specific height, when the droplet has reached the lower threshold contact angle θ min the contact line depins.This causes a distinct corner in the H(t) curve as the linear relation between height and volume gives place to the usual cubic relation V ∼ R 3 in the diffusion-limited case (Fig. 6 (a)) and a linear relation with a different slope in the transition-limited case (Fig. 6 (c)).
The constant evaporation rate in the case of a pinned contact line also appears in the approximate expressions of the total evaporation rate derived in appendix A. There, the evaporation rates for both, the strongly diffusion-limited and strongly transition-limited case, depend only on the radius R, but not on any other characteristics of the drop shape.
In the transition-limited case, the evaporation rate is quite uniform over the drop surface (∼ R 2 in the long-wave limit), whereas in the diffusion limited case, evaporation mostly occurs in the vicinity of the contact line, i.e., it scales with the droplet radius R.
Once the contact line has depinned, the droplet continues to shrink with a contact angle close to θ min , where the curves for volume V (t) and height H(t) recover the characteristic shape of the case without contact angle hysteresis studied above.Note, in particular, that their slope matches exactly for the transition-limited case in Fig. 6 (c).
Overall, due to the constant rate of evaporation, in all cases pinning results in a faster drop evaporation.The effect is, however, most pronounced in the phase transition-limited case, where the evaporation rate exhibits a dependency on the drop surface area ∼ R 2 .The effects are more subtle for the strongly diffusion-limited case, where the drop shape is mostly irrelevant as the vapour diffusion is the limiting factor.

Comparison to Stokes description and experiments
Next, we compare the presented results obtained with the developed thin-film model to simulations with the Stokes model in a gap geometry and to experiments with evaporating sessile drops in a narrow gap.In Fig. 7, snapshots obtained with the Stokes model are shown, at parameters and times that exactly coincide with the ones employed in Fig. 3 for the thin-film model.In the diffusion-limited case, Fig. 7 (a), the evaporation rate is highest near the contact line.However, in contrast to the thin-film model, where the Derjaguin pressure g (h) has an influence on the evaporation rate, this corresponding effect is reflected in the Stokes model by the change from drop to bare substrate.With other words the jump from a finite evaporation rate at the drop edge to zero at the bare substrate is slightly smoothed by the Derjaguin pressure in the thin-film model.A stronger difference can be seen in the vapour distribution near the droplet apex that here clearly shows a deviation from the z-independent distribution assumed in the thinfilm model.However, a distance away from the droplet, the vapour distribution indeed becomes more and more uniform in the direction vertical to the substrate, i.e., diffusion turns into a simple one-dimensional radial process as also encoded in the thin-film model.The rate-limited case in Fig. 7 (c) features a uniform evaporation rate and a vapour diffusion that is sufficiently fast to approach an almost homogeneous vapour concentration fixed by the lab vapour concentration.This is in full agreement with the corresponding results of the thin-film model in Fig. 3 (c).
Next we assess the influence of the approximations made in the mesoscopic model in the shallow-drop case, i.e., assuming a z-averaged vapour distribution, a simplified expression for the curvature resulting in parabolic drop profiles, and laminar flow.To do so Fig. 8 gives the volume and height evolution for both models.Note that, due to the parabolic drop profile in the thin-film model and the spherical-cap shaped profile in the Stokes model, the equilibrium contact angle of 44.38°(as used for the thin-film results) had to be slightly reduced to 43.38°in the Stokes description to simultaneously have identical initial heights and volumes in the two models.
Obviously, for the diffusion-limited case in Fig. 8 (a), the agreement is better for smaller gap heights, i.e., for smaller distances between the droplet and the upper plate.For d = 1.1 mm, results in the initial stage of the evolution almost coincide, whereas they deviate towards the final stage, when the droplet height h becomes considerably smaller than the plate distance d.The diffusive vapour transport in the Stokes model is slower than predicted by the purely radial transport assumed in the thin-film model.The larger the initial distance of drop apex to upper plate, the higher the initial deviation.For the largest considered plate distance of 10 mm, the assumption of a z-averaged vapour distribution underpredicts the total evaporation time by a factor of 3.
In the intermediate case, Fig. 8 (b), it is apparent that at small gap widths the Stokes model shows nearly identical but slightly faster evaporation than the thin-film model.Since such slightly faster evaporation is also visible for the transition-limited case in Fig. 8 (c), it can be attributed to the different mathematical treatment of the free surface: in the shallow-drop case, the evaporation rate is effectively determined by the base area πR 2 of the parabolic drop.In the Stokes model, the full surface area of a spherical drop is considered, leading to an increased total volume loss rate once the transfer rate becomes the limiting factor.Note that this discrepancy disappears in the full-curvature formulation of the mesoscopic model, implying improved expressions for Laplace pressure and evaporation terms.For large plate distances d in the intermediate case, however, our central approximation of a z-averaged vapour density again results in an underprediction of the total evaporation time in the mesoscopic model.
The above argument is further confirmed by the transition-limited case in Fig. 8 (c).There, the macroscopic Stokes calculations uniformly give a slightly faster evaporation than in the mesoscopic model for all considered gap heights.This indeed implies that the cause of the difference is in the treatment of the drop surface area as discussed before.In the full-curvature formulation, mesoscopic and macroscopic approaches fully coincide.
We conclude that as expected the agreement of the results obtained with the developed mesoscopic approach and the macroscopic Stokes approach improves for smaller gap heights d, as then the gas phase becomes increasingly vertically homogeneous also in the Stokes simulations.Similarly, smaller contact angles θ e allow for lower drop heights and thus for lower gap width, therefore improving the agreement.
Finally, we compare the results obtained from the Stokes and thin-film simulations to the experiments.Therefore we adapt the parameters and initial condition of the simulations such that they closely match the experiment.The resulting parameters are the same as mentioned in section 2.2.4,with the only difference that we here assume a relative humidity of ρ lab /ρ sat = 50 % at the far end of the gap (at r = 25 mm).Additionally, in the simulations the phase transition coefficient is chosen as M = 4 × 10 −10 m/Pa s, i.e., the process is limited by diffusion rather than phase change.The experimental drop shape and volume were extracted from the shadowgraphy images and used to initialise the simulations with droplets of identical size (see table 1 in the appendix for the measured initial drop shape parameters).Since in the experiment the droplet exhibits a strong pinning effect at its initial contact line position and depins only at a very low contact angle (see also Fig. 2), we here pin the contact line in the thin-film model using the mechanism described in section 3.2 and replace the Navier slip condition with a no-slip condition in the Stokes model.
In Fig. 9 we depict the normalised volume and height curves as obtained from experiments and simulations for varying values of the gap height d.The gap height is color coded, while the line type distinguishes the method (thin-film model / Stokes model / experiment) used for the generation of the data.It should be noted that the case of having no upper plate at all (d → ∞) is included only for the experiment and Stokes simulations as this limit cannot be achieved in the thin-film model.
As discussed previously, the thin-film and Stokes model agree well, in particular for small gap heights, and their discrepancies are explained by the effect of vertical diffusion in the Stokes model.Further, we find that both models compare well to the experimental results, where again the agreement is best for smaller heights d.Both, in the experiment and in the simulations, the evaporation rate is higher for larger values of the gap height and the drop evaporates most rapidly in the unconfined case d → ∞.This reveals that the experiment is in fact in a diffusion-limited regime, as are the simulations.The volume is found to decrease linearly over time both in the experiment and in the modelling, as it is expected for the diffusion-limited case.It is again due to the contact line pinning that the drop height H relates linearly to the volume, i.e., H(t) is also linear.However, this relation relaxes in the experimental data once the drop becomes very shallow and the contact line depins.The depinning causes a deformation of the H(t) curve that Comparison between experiments (crosses), Stokes simulations (dashed lines), and the thin-film model (solid lines) in the diffusion limit for varying gap height d, or no gap at all (d → ∞, unfeasible in the thin-film model).The modelling parameters were adapted to the experiment, their values are noted in the main text.Here, the contact line is pinned in both simulations and experiments but is found to depin in some experimental runs.The agreement between Stokes and thin-film simulations is best for small gap heights, as discussed in the context of Fig. 8. Compared to the experiment both simulations underpredict the evaporation time, presumably due to the exclusion of thermal effects in the models.
resembles exactly our simulations where contact angle hysteresis was included in the thin-film model, namely, Fig. 6.
Despite the similarities between experiment and theory, the modelling appears to systematically overestimate the evaporation rate even before the depinning occurs.We presume that this is caused by a reduction of the vapour pressure due to cooling by latent heat of evaporation.As both presented models are isothermal, such effects are excluded in the simulations.This is supported by Figure 12 in the appendix, where we showcase the same measurements performed with an extended Stokes model that includes heat transfer [model not discussed here, specification follows Diddens (2017) for the case of a simple liquid].There, most of the discrepancies are explained by the inclusion of thermal effects.The assumption of an isothermal gas phase is more accurate for thinner gap heights, hence the agreement between the isothermal models and experiments (or thermal modelling) is best for small gaps.

Conclusion
We have developed an effective two-field thin-film model for sessile shallow drops of volatile liquids under the strong confinement imposed by placement in the narrow gap between two parallel plates.The model consists of coupled time evolution equations for the drop height or film thickness profile and the vertically averaged vapour concentration profile in the narrow gap.A particular strength lies in its ability to describe the coupled liquid convection, liquid-vapour phase change and vapour diffusion for the full spectrum of dynamical coefficients ranging from the diffusion-limited regime at fast phase change to the phase transition-limited regime at slow phase change.In this way it unifies the two groups of isothermal thin-film models available in the literature as reviewed in the introduction.
After presenting the model it has been used to investigate evaporating drops in the different regimes.The obtained numerical results have on the one hand recovered analytical relations we have derived in the two limiting cases.On the other hand they have provided insights into the intermediate regime where the evolution shows aspects of both limiting cases.Furthermore, additionally to the thin-film model we have introduced a model based on the Stokes equation for the liquid in the drop, a diffusion equation for the vapour concentration in the gap and adequate boundary conditions.In particular, the employed modelling of the process of phase transition matches the one in the thinfilm model.Results obtained with the thin-film model and with the Stokes model are compared and, overall, show satisfactory agreement.As expected, the agreement is nearly perfect in all dynamic regimes for narrow gaps and relatively large drops and becomes less so for wide gaps and small drops.
We found that in the long-wave model the dependency on the contact angle can be traced back to a dependency on the liquid surface area.Thereby, the evaporation rate is independent on the contact angle in the diffusion-limited case whereas in the phase transition-limited case the evaporation of drops with smaller contact angle (or larger initial radius R) speeds up by 1/R 2 .Note that the former is only true in the constrained gap geometry and not in an unconfined setting.Additionally, the reduced height of smaller contact angle droplets allows for more shallow gap geometries, which improves the agreement of the long-wave results with the full Stokes model.
The good agreement with the Stokes results for narrow gaps for small or moderate contact angles indicates that the developed thin-film model can be used as a valid tool for the study of shallow sessile drops of volatile liquids in a variety of contexts across all dynamical regimes.The model can easily be adapted for many closely related systems, e.g., incorporating chemically or topographically heterogeneous substrates, resulting in contact line pinning and depinning (van Der Heijden et al. 2018), adding gravity or other body forces (Kumar et al. 2020), or considering the importance of the evaporation regime in dip coating (Bindini et al. 2017;Dey et al. 2016).As the presented thin-film description is of gradient dynamics form, necessary changes can often be incorporated via adapting the underlying free energy functional that here only incorporates wetting, interface and bulk energies of the liquid and vapour entropy [cf.Eq. (2.9)].This includes the usage of other equations of state (then adapting the considerations in Section 2.2.2), employing the full (not long-wave) curvature for the liquid-vapour interface (Thiele 2018), and the incorporation of more realistic wetting energies (Churaev 2003;Hughes et al. 2017).
Note that the presented thin-film model is isothermal, while the Stokes calculations are presented for both cases, isothermal and non-isothermal.One may incorporate thermal effects into the thin-film model, in the simplest case by using an effective film height-dependent transfer constant as often done for drops on heated substrates in a gas of low thermal conductivity [see discussion of Eq. ( 11) in (Thiele 2014)].Additionally, including thermal Marangoni effects as in (Ajaev et al. 2010;Savva et al. 2017) and also heat diffusion is challenging as the full dynamics of heat has to be incorporated into the gradient dynamics form.This remains a task for the future.However, because of its gradient dynamics form, the presented thin-film model may readily serve as a building block in thin-film descriptions of more complex systems.It may, for instance, be combined with recent descriptions of droplets on soft elastic substrates (Henkel et al. 2021(Henkel et al. , 2022) ) and polymer brushes (Thiele & Hartmann 2020) to study how the evaporation or condensation of drops is modified by substrate softness and brush properties, respectively.In principle, the approach can also be extended beyond pure liquids to capture evaporation, condensation and absorption of liquid mixtures, e.g., to study the dynamics of mixtures of volatile liquids where selective evaporation results in intriguing effects (Christy et al. 2011;Karpitschka et al. 2017;Hack et al. 2021), the evaporation of drops on liquid-infused (Guan et al. 2015) or porous (Gambaryan-     1. Characteristic parameters of the initial drop shape (t = 0) in the experimental runs for the comparison with simulations in Fig. 9.The values are measured from shadowgraphy images and are then used to initialise the simulations with matching drops.

SFigure 1 .
Figure1.Sketch of the considered system.It consists of two parallel rigid smooth solid plates separated by a gap of width d.A layer or shallow sessile drop of volatile liquid with a thickness profile h(r, t) with r = (x, y) T is situated on the lower plate without touching the upper one.The height-averaged density between the liquid and the upper plate is ρ(r, t) for the vapour, and ρair(r, t) for the remaining (inert) gas, i.e., the total height-averaged gas density is ρtot = ρ+ρair.There can be exchange between the liquid and the vapour due to evaporation/condensation.
.31) These conditions also ensure zero liquid and vapour flux through the computational domain boundary at r = 0.At the outer boundary far from the drop the gap between the plates is open, i.e., air and vapour can freely be exchanged with the surrounding laboratory.Thus, we assume a constant lab humidity (or vapour concentration) ρ lab corresponding to a Dirichlet BC at r = L. Together with natural Neumann BC for the film profile h and the chemical potential δF/δψ 1 , this gives the remaining BC ∂ r h r=L = 0, ∂ r δF δψ 1 r=L = 0, and ρ r=L = ρ lab .

Figure 2 .
Figure 2.An example of experiments of sessile droplet evaporation in a gap.The red solid line is the bottom plate.The blue dashed lines shows the top plate.The initial shape of the droplet is shown with a dotted line in all panels.For this test, the initial diameter of the drop was 1.33 mm and the gap size was d = 0.88 mm.

Figure 3 .
Figure 3. Single simulation snapshots are shown obtained with the thin-film model (2.30).Panel (a) gives the diffusion-limited case, (b) an intermediate case, and (c) the phase transition-limited case at evaporation rate constant M as given above the panels.The gap height is fixed at d = 3 mm.All other parameters are identical and correspond to the values given at the end of Section 2.2.4.The top panels show the local evaporation rate jevap(r, t) while the bottom panels give the droplet profile (solid line, filled blue, left hand side scale) and the local relative humidity (dotted line and bluish shading, right hand side scale).The inset in the top panel (a) magnifies the peak in evaporation rate at the contact line (lateral size 0.7 mm).

Figure 4 .
Figure 4.The thin-film model is used to determine the change in (top) drop volume and (bottom) drop height for various gap heights d given in the legend in (a) the diffusion-limited case, (b) an intermediate case, and (c) the phase transition-limited case via different values of the phase transition coefficient M .The inset in (b) shows the same data but shifted in time such that the curves coincide at the time the drop has fully evaporated, revealing that the curves show the same phase transition-limited behaviour when the droplet is much smaller than the gap height d.

Figure 5 .
Figure5.Droplet volume and height over time for long-wave simulations identical to Fig.4but with varying values of the contact angle parameter θe.The gap height is now fixed at d = 1.5 mm.All drops were initialised with their respective contact angle, same radius but varying initial volume V0.The normalised time scale t/V0 accounts for the different evaporation times of the different drop volumes, revealing that the contact angle itself has little to no qualitative effect on drops of the same radius.

Figure 6 .
Figure6.Droplet volume and height over time for long-wave simulations identical to Fig.4but with an initial pinning of the drop's contact line due to a contact angle hysteresis of 30°.For comparison, the simulations without pinning are included as thin dashed lines.The noticeable corner in the height curves indicate the depinning event.While the contact line is pinned, the volume evaporation rate dV /dt remains nearly constant, accelerating evaporation particularly in the transition-limited case.

MFigure 7 .
Figure 7. Simulation snapshots obtained using the Stokes model at identical parameters and times as employed in Fig. 3 for the thin-film model.Panel (a) gives the diffusion-limited case, (b) an intermediate case, and (c) the phase transition-limited case.The evaporation rate is indicated by the arrows at the liquid-vapour interface, whereas the vapour density ρ(r, z, t) is given by the background shading.
The intermediate case depicted in Fig. 7 (b) shows similar features as the corresponding thin-film model result in Fig. 3 (b), i.e. an almost uniform evaporation rate with a slight enhancement near the contact line, where the vapour diffusion still contributes as limiting factor, i.e. vapour gradients are still visible.S. Hartmann, C. Diddens, M. Jalaal and U. Thiele = 4 × 10 −10 m/Pa s b) M = 4 × 10 −16 m/Pa s c) M = 4 × 10 −18 m/Pa s

Figure 8 .
Figure 8.Comparison of the volume evolution V (t) and the height evolution h(t) for the thin-film model (solid lines) and the Stokes model (dashed lines) for (a) the diffusion-limited case, (b) an intermediate case, and (c) the phase transition-limited case.
Figure9.Comparison between experiments (crosses), Stokes simulations (dashed lines), and the thin-film model (solid lines) in the diffusion limit for varying gap height d, or no gap at all (d → ∞, unfeasible in the thin-film model).The modelling parameters were adapted to the experiment, their values are noted in the main text.Here, the contact line is pinned in both simulations and experiments but is found to depin in some experimental runs.The agreement between Stokes and thin-film simulations is best for small gap heights, as discussed in the context of Fig.8.Compared to the experiment both simulations underpredict the evaporation time, presumably due to the exclusion of thermal effects in the models.

Figure 11 .
Figure 11.Highly zoomed snapshot of the contact line region of Fig. 3 (a) in order both showing the character of the peak in the local evaporation rate hevap and the smoothness of the solution in the long-wave model.The apogee of the evaporation rate smoothens over a very short length, as seen in the inset.

Figure 12 .
Figure12.Comparison of the isothermal Stokes model used in this work (dashed lines) and an extended Stokes model that includes heat transfer (solid lines) to the experimental data (crosses), proving that the thermal effects can explain most of the discrepancies between our isothermal models and the experiments. gap height d drop height drop volume contact angle drop radius