Does the van der Waals force play a part in evaporation?

It is argued that the van der Waals force exerted by the liquid and vapor/air on the molecules escaping from one phase into the other strongly affects the characteristics of evaporation. This is shown using two distinct descriptions of the van der Waals force: the Vlasov and diffuse-interface models, each of which is applied to two distinct settings: a liquid evaporating into its vapor and a liquid evaporating into air (in all cases, the vapor-to-liquid density ratio is small). For the former setting, the results are consistent with the Hertz--Knudsen Law (HKL), but the evaporation/condensation probability is very small (in the classical HKL, it is order one). For the latter setting, the dependence of the evaporation rate on the difference between the saturated vapor pressure and its actual value is shown to be nonlinear (whereas the classical HKL predicts a linear dependence). The difference between the two settings indicates that the van der Waals force exerted by the air strongly affects evaporation (contrary to the general assumption that the ambient gas is unimportant). Finally, the diffuse interface model is shown to be inapplicable in a narrow region at the outskirts of the interface -- as a result, it noticeably underestimates the evaporative flux by comparison with the (more accurate) Vlasov model.


I. INTRODUCTION
Evaporation is fundamental to numerous environmental, biological, and industrial processes, and the Hertz-Knudsen Law (HKL) is our primary tool for modeling it.This paper argues that the HKL needs to be modified by taking into account the effect of the van der Waals force on the evaporative flux.
In its original formulation 1,2 , the HKL was based on an assumption that the flux of molecules escaping from a liquid into vapor does not depend on the vapor pressure -hence, this flux can be calculated as if the vapor were saturated.Calculating also the flux in the opposite direction (which does depend on the actual vapor pressure), one can show that the net evaporative flux is where ρ (v.sat) is the saturated vapor density, ρ (v) the actual density, R the specific gas constant, and T the temperature (assumed, for simplicity, to be the same in the liquid and vapor).
Expression (1) does not involve a single adjustable parameter and, thus, is unlikely to be accurate for all liquids under all conditions.To make it more adaptable, one can assume that some of the escaping molecules bounce back, as do those travelling in the opposite direction.It can be argued 3 that a molecule's probability of evaporation equals that of condensation, resulting in the following modification of expression (1): a) Email address: Eugene.Benilov@ul.ie;Homepage: https://eugene.benilov.com/where the evaporation/condensation probability θ (called also "mass adjustment coefficient") depends on T .The amended version of the HKL still disagrees with some of the available experiments, and those disagree with each other: for, say, water, the measured values of θ vary between 0.01 and 1 for the same temperature 4,5 .There are also several theoretical models (e.g., [6][7][8] ), but the discord in the experimental results makes it difficult to choose the most accurate theory.The present paper is motivated by an observation that none of the existing models of evaporation accounts for the van der Waals (vdW) force.Yet it is clearly important: it holds the liquid/vapor interface together (by balancing the pressure gradient due to the density variation) -hence, should affect the molecules passing through the interface.
It can also be argued that the vdW force makes evaporation of a liquid into its vapor different from evaporation into air.To understand why, note that the vdW force exerted by the bulk of the liquid pulls the escaping molecules back and, thus, impedes evaporation -while the vapor and air pull them forward and, thus, encourage evaporation.Since under normal conditions the vapor and air densities differ by orders of magnitude, the former exerts a much stronger vdW force than the latter.As a result, evaporation into air should occur much faster than that into vapor -and this is one of the two main conclusions of the present work.
The other one is less intuitive, but still has important physical implications.As shown for evaporation into air, the vdW force makes the dependence of the flux E on the density difference ρ (v.sat) − ρ (v) nonlinear, whereas the HKL predicts that E ∼ ρ (v.sat) − ρ (v) [see Eq. (2)].The two results can be reconciled only if the evaporation/condensation probability θ in (2) depends on ρ (v) .Such a dependence could explain the above-mentioned discord in the measurements of θ for the same liquid at the same temperature.
The present paper employs two different descriptions arXiv:2312.11117v2[cond-mat.soft]12 Feb 2024 of the vdW force: the Vlasov model and the diffuseinterface model.The former has been used before to study interfaces 9 , contact lines, and liquid films [10][11][12] -but not evaporation; the latter has been applied to evaporation [13][14][15] , but its connection to the HKL has not been properly explored.
The two models of the vdW force will be used in conjunction with the isothermal Navier-Stokes equations.This simple framework is sufficient for demonstrating the importance of long-range intermolecular forces for evaporation in principle, and this is the aim of the present paper.
In what follows, Secs.II-III examine evaporation of a liquid into its vapor, Secs.IV-V examine evaporation into air, and Sec.VI explains why the (alleged) shortcomings of the HKL have not been so far observed by the experimentalists and researchers working on molecular dynamics.
Since the material presented in this paper is associated with a number of bulky specialized terms, several abbreviations will be used.For future reference, they are listed in Table I.

II. EVAPORATION OF A LIQUID INTO ITS VAPOR: THE FORMULATION A. Thermodynamics
A model of phase transitions should account for the fluid's thermodynamic properties.These are described in this subsection, in a brief but self-consistent manner.
Following 16 , one can fully characterize a fluid by setting the dependence of its specific (per unit mass) internal energy e and entropy s on the density ρ and temperature T .The functions e(ρ, T ) and s(ρ, T ) are not arbitrary, but are constrained by the so-called Gibbs fundamental relation, which can be written in the form (the equivalence of this equality to the traditional form of the Gibbs relation is demonstrated in Appendix A of Ref. 17 ).Given e(ρ, T ) and s(ρ, T ), one can find the pressure p(ρ, T ) and chemical potential G(ρ, T ) via the formulae Eqs. ( 4)- (5) imply that which is the only thermodynamic identity needed in this paper.(10) for water at T = 352 • C. The region where the vapor density does not have a match in the liquid region is shaded (it exists only if T is sufficiently high, so that the local minimum of the curve p vs. ρ is above the horizontal axis).
An example of p(ρ, T ) (often referred to as the equation of state) is shown in Fig. 1.Observe that the dependence of p on ρ is nonmonotonic: the states between the origin and local maximum are vapor and those between the local minimum and infinity, liquid.The states between the minimum and maximum (those with ∂p/∂ρ < 0) are unstable.To understand why, observe that, in this case, a spatially-localized density increase causes a pressure decrease -the resulting pressure gradient generates an inward flow -which causes a further density increasehence, instability.
The saturated vapor density ρ (v.sat) and the matching liquid density ρ (l.sat) , for a given temperature T , satisfy the so-called Maxwell construction -i.e., p(ρ (v.sat) , T ) = p(ρ (l.sat) , T ), ( 7) Conditions ( 7)-( 8) guarantee that the interface separating the vapor and liquid is in mechanical and thermodynamic equilibrium, respectively.Due to nonmonotonicity of p, the solution of Eqs. ( 7)-( 8) is not unique.To make it such, require that the vapor and liquid are both stable The general results reported in this paper are illustrated using the so-called Enskog-Vlasov fluid model, according to which where c is the specific heat capacity of the fluid under consideration.
The first term in e and the first two in s correspond to ideal gas.The term involving a describes the contribution of the vdW force to the internal energy, and the function Θ(ρ) is the non-ideal contribution to the entropy.Both a and Θ(ρ) should be fixed by fitting the fluid's equation of state to its empiric shape 17 ; in the present paper, the values corresponding to water will be used (see Appendix A 1).The Enskog-Vlasov fluid model is sufficiently accurate (as shown in Secs.8.1-2 of Ref. 17 ), plus it is consistent with the Vlasov description of the van der Waals force used in this paper.
Substituting the Enskog-Vlasov expressions (9) into Eqs.( 4)-( 5), one obtains where • • • hides the terms in the chemical potential depending only on T (they cancel out from identity (6), the Maxwell construction, and all the equations to come).In the low-density limit, expressions ( 10)-( 11) yield These asymptotics describe an ideal gas and, thus, are not specific to the Enskog-Vlasov fluid model.The first term in expression (10) will be referred to as the thermal pressure; denoting it by p, one obtains p = p + aρ 2 . ( A similar equality inter-relates the thermal and full chemical potentials, Ĝ = G + 2aρ. Eqs. ( 13)-( 14) and ( 6) imply that Finally, the low-density asymptotics of the thermal chemical potential coincides with that of the full chemical potential, Ĝ ∼ RT ln ρ as ρ → 0, (16)   whereas the asymptotic of p will not be needed.

B. The low-temperature limit
Assume that the nondimensional temperature T nd is small, This restriction typically applies to all common liquids under normal conditions; for water between 0 • C and 100 • C, for example, one obtains 0.065 ≲ T nd ≲ 0.082.
The assumption T nd ≪ 1 implies that the vapor-to-liquid density ratio is also small -or, to be precise, the smallness of T nd makes ρ (v.sat) /ρ (l.sat) exponentially small (as shown in Ref. 18 for the generic van der Waals equation of state, but is also true generally).For water between 0 • C and 100 • C, for example, one can use the online calculator 19 to obtain 4.9 × 10 −6 ≲ ρ (v.sat) ρ (l.sat) ≲ 6.2 × 10 −4 .

C. The van der Waals force
Consider a compressible fluid characterized by its density field ρ(r, t), where r is the position vector and t, the time.Introduce also the molecular mass m, so that ρ/m is the number density.
Let the potential of the van de Waals force exerted at a point r by a molecule located at a point r ′ be m 2 Φ(r−r ′ ), where the factor m 2 is introduced for convenience.If the fluid is isotropic, then Φ(r) = Φ(r), where r = |r|.Without loss of generality, one can assume that Φ → 0 as r → ∞.
Strictly speaking, Φ(r) should be determined by examining quantum interaction of the fluid's molecules.In practice, however, microscopic characteristics like Φ(r) are determined by calculating the corresponding macroscopic parameters of the fluid and fitting them to their empiric values.Let the fluid occupy the whole space, so that the volumetric density of the collective force, induced by all the molecules at a point r, is where the integral in (18) is to be evaluated over the whole space (the same is implied in all further integrals with omitted limits).Expression (18) will be referred to as the Vlasov model, which is how the collective field approach is called in plasma physics 20 .
If ρ depends only on the vertical coordinate z (i.e., physically, the liquid/vapor interface is flat and horizontal), the integral in ( 18) can be rewritten in cylindrical coordinates (r ⊥ , β, z ′ ) and reduced to where The diffuse-interface model (DIM) is based on the assumption that the spatial scale of Ψ(z) [inherited from the original intermolecular potential Φ(r)] is much shorter than that of the density field.Then, the integral on the right-hand side of ( 19) can be simplified by changing z ′ → z − z ′ , expanding ρ(z − z ′ , t) in powers of z ′ , and truncating the expansion after the first three terms.The integral involving z ′ Ψ(z ′ ) vanishes (because Ψ(z ′ ) = Ψ(−z ′ ) due to isotropy), and one obtains 21,22 where the factor 2 is introduced for convenience, and will be referred to as the van der Waals parameter and Korteweg parameter, respectively.The former is the same as its namesake a in the Enskog-Vlasov fluid model ( 9)-( 11) -hence, its value for a particular fluid can be deduced by fitting the Enskog-Vlasov equation of state to its empiric counterpart (for water, this is done in Appendix A 1).The value of the Korteweg parameter K, in turn, can be deduced from the fluid's capillary properties (for water, see Appendix A 2).
Even though a and K were introduced as parameters of the DIM, they can be viewed as global characteristics of the general Vlasov potential Φ(r).Furthermore, since they are integral characteristics of Φ(r), they are more important than its actual shape.In particular, Eqs. ( 22)- (23) suggest that the spatial scale of the vdW force is which is one of the crucial characteristics in interfacial dynamics.
In the present paper, the following example of the Vlasov potential is used where H(∆ − r) is the Heaviside step function, and B, C, and ∆ are adjustable constants.Substituting ( 25) into (20), one obtains To adapt Ψ(z) to the fluid under consideration, one should express B and C through the van der Waals and Korteweg parameters, a and K. Substituting expression (26) into ( 22)-( 23) and solving for B and C, one obtains With a and K deduced from empiric data, the undetermined parameter ∆ can be viewed as the one describing the shape of Φ(r).In particular, Φ(r) is monotonic only if If ∆ is outside this range, the vdW force (which is ∼ ∇Φ) is repulsive for some r -whereas physically, it should be attractive.Thus, ∆ should better be chosen from range (29).
One might also argue that the vdW force should not involve a small-r component: the short-range part of the intermolecular interaction is responsible for collisions and supposed to be accounted for by the viscosity term (in hydrodynamics) or collision integral (in kinetic theory) -not the Vlasov term.With this in mind, one should make ∇Φ near r = 0 as small as possible; in terms of expression (25), this corresponds to B = 0, so that ( 27) yields Evidently, this value is included in range (29) as its lefthand endpoint.
Since particular case (30) satisfies all the criteria, it will be used in the remainder of this paper.
The following comments are in order: • Even if one chooses a different ∆ from range (29), the fluid's macroscopic properties remain virtually the same.In particular, when a and K are fixed while ∆ varies through the whole range (29), the corresponding change of the surface tension γ is approximately 0.1% (this calculation was carried out using the Enskog-Vlasov equation of state with the parameters of water at 25 • C; see also Appendix A 2 for the dependence of γ on the equation of state, a, and K).
• The choice of the formula for Φ(r) appears to also be unimportant: an exponential alternative to (25)  was tested, and the dependence of γ on the temperature was found to be virtually the same.
Note that the diffuse-interface model corresponds to the limit ∆ → 0 and, thus, is not included in interval (29).Another problem with the derivation of the DIM via expansion (21) has been noted in Ref. 11 : for some intermolecular potentials, the higher-order terms omitted from (21) involve divergent integrals.For potential (25), this problem does not arise -but the mere possibility of a divergence may indicate "a qualitative difference between the solutions of the exact and truncated equations" 11 .
All this does not mean that the DIM should be discarded; it has been used for more than a century by hundreds of researchers -and deserves to be examined at face value.This is what is done in the present paper, and the results obtained are tested against those of the more general Vlasov model.

D. Governing equations
Evaporation of common liquids at normal conditions is a slow process -hence, the slow-flow approximation can be used, which amounts to the following equations: where w(z, t) is the vertical velocity, p is the thermal pressure, and µ is the effective viscosity related to the shear and bulk viscosities by The non-thermal part of the pressure comes from the van der Waals force: in the DIM, this can be shown explicitly (by substituting expression (21) for F into the momentum equation (32) and incorporating the term ∼ a into the pressure term).In the general Vlasov model, however, the non-thermal pressure has to remain 'hidden' inside the force term.The viscosity µ(ρ, T ) should be treated as a known function.
4][15] for the DIM, evaporation of common liquids under normal conditions is not sensitive to the finite-ρ range of the viscosity function µ(ρ, T ); its only important characteristic is the low -density limit µ 0 .The same is true for the Vlasov model (more details given below) -thus, when illustrating the general results, the simplest (density-independent) approximation of the viscosity will be used, µ(ρ, T ) = µ 0 (T ).For specific computations, one still needs the dependence of µ 0 on T -in this paper, that of water is used (see Appendix A 3).

E. Boundary conditions
Consider a flat horizontal interface separating a liquid and its vapor (the former is located below the latter).Mathematically, this corresponds to where ρ (v) and ρ (l) are the vapor and liquid densities, respectively, and the z axis is directed upwards.Note that ρ (v) is a given parameter (determined by the relative humidity), whereas ρ (l) is to be found.If the vapor is undersaturated (ρ (v) < ρ (v.sat) ), evaporation gives rise to a flow -such that i.e., physically, the liquid far below the interface is at rest, and the vapor flow far above the interface is stress-free.Equations ( 19), ( 31)- (32) and conditions (34)-(37) form a boundary-value problem to be solved.

F. The isobaricity condition
Under the DIM, one can readily show that equation (32) and boundary conditions (34)- (37) imply that the pressure values far above and far below the interface are equal, i.e., p(ρ (l) , T ) = p(ρ (v) , T ).(38)  For the Vlasov model, this result holds too, but is a little harder to prove (see Appendix B).Equality (38) will be referred to as the isobaricity condition; it allows one to calculate the density ρ (l) of the evaporating liquid without solving the governing equations.To do so, recall that the vapor density ρ (v) is a known parameter, and treat (38) as an equation for ρ (l) .
The isobaricity condition has another important implication: if the temperature is sufficiently high and/or the vapor density ρ (v) is sufficiently low, Eq. ( 38) does not have a solution for ρ (l) -see an illustration in Fig. 1.In such cases, the liquid below the interface cannot be homogeneous; it was conjectured in Ref. 14 that it boils, but perhaps this phenomenon should be called cavitation.
Either way, this effect will not be discussed in further detail (because it typically occurs at a temperature much higher than "normal").One should only remember that solutions describing steady evaporation may cease to exist when the temperature exceeds a certain threshold.

III. EVAPORATION OF A LIQUID INTO ITS VAPOR: THE SOLUTION
Assume that evaporation is steady -hence, the liquid/vapor interface recedes at a constant velocity equal to −E/ρ (l) , where E is the evaporation rate and ρ (l) , the liquid density.Thus, seek a solution of the form where In terms of the new variable, the density equation ( 31) becomes (the subscript new omitted) Integrating this equation and fixing the constant via boundary conditions ( 34) and ( 36), one obtains Coincidently, this expression satisfies the second boundary condition for w, (37).
Next, rewrite the momentum equation (32) in terms of z new , omit new , use (39) and (19) to eliminate w and F , then use identity (15) to express p through Ĝ, and eventually obtain This equation and boundary condition ( 34)-( 35) are invariant with respect to an arbitrary shift z → z + const -hence, they do not fully fix the solution ρ(z).An extra boundary condition is needed -say, where the saturated densities are used to 'pin' solutions with different ρ (v) to the same point of space associated with the equilibrium solution.

A. The diffuse-interface model
The DIM equation for steady evaporation can be obtained similarly to its Vlasov counterpart, Eq. ( 40) -one only needs to represent the vdW force by the differential expression ( 21) instead of its integral counterpart (19).Eventually, one obtains Boundary-value problem (34)-( 35), ( 41), ( 43) was solved numerically using the function BVP5C of MATLAB.A typical solution is shown, together with the equilibrium solution, in the upper panel (labeled "DIM") of Fig. 2.
The following features should be observed: • Within the interface, the equilibrium and nonequilibrium solutions are indistinguishable.
• The two curves split when passing through a narrow region just outside the interface (to be referred to as the van der Waals layer ) and remain constant after that.
These observations help one to examine the problem asymptotically in the low-temperature limit.Two asymptotic zones can be identified: the interface and vdW layer.In the former, the solution is determined by the balance of the pressure gradient and van der Waals force, whereas in the latter, viscosity comes into play.Since the interface is, essentially, in equilibrium, it is the vdW layer that determines the evaporation rate.
The asymptotic solution of boundary-value problem (34)-( 35), ( 41), ( 43) is described in Appendix D, and is summarized here in terms of the temperature T and relative humidity The equilibrium interface (marked with "=") and nonequilibrium interface with H = 0.5 (marked with "∼"); in both cases, T = 50 • C. The results in panels labeled by "DIM" and "VM" are computed using the diffuse-interface and Vlasov models, respectively, with the parameters of water (also implied in all further figures).
As shown in Appendix D, the evaporation rate is where is of the same dimension as E, whereas ẼD (H) is nondimensional [and determined by boundary-value problem (D15)-(D17)].
The exact and asymptotic solutions, both found numerically, are illustrated in the upper panels (labeled "DIM") of Fig. 3.The following features should be observed: • Panel DIM(a) shows that the asymptotic solution becomes inapplicable near the point where the exact solution ceases to exist.This comes as no surprise, as this temperature is fairly close to the critical point -hence, ρ (v.sat) /ρ (l.sat) is not small there.
• Panel DIM(b) shows that the exact curves for different T collapse, or nearly collapse, onto the same asymptotic curve.This is a result of the 'separation' of T and H in expression (44).
• The solid curves in panel DIM(b) are not extended to small H because the exact solution is difficult to compute in this region.The difficulty is caused by the smallness of the vapor density, so that the last term in Eq. ( 43) becomes nearly singular.
• The asymptotic curve in panel DIM(b) is close to a straight line (although it is not one exactly).
As shown in Appendix D, the spatial scale of the vdW layer, l L , and that of the vdW force, l F [given by (24)], are such that where T nd is the nondimensional temperature defined by (17).Recalling that, for common liquids under normal conditions, T nd is small, whereas ρ (v.sat) /ρ (l.sat) is exponentially small, one concludes that l L ≪ l F .This is the opposite of what the DIM is applicable to.This observation motivated the author of the present paper to re-examine the problem using the Vlasov model (whose applicability does not require that the flow's spatial scale be large).

B. The Vlasov model
The Vlasov equation ( 40) is much harder to solve numerically than its DIM counterpart (43): the former is of an integro-differential kind, for which there are no readymade tools in MATLAB or similar packages.The only numerical algorithm the author of this paper has come up with (see Appendix C) does not perform well for small H and/or T , and often requires manual fine-tuning of the computational parameters.
A typical solution of Eq. ( 40) subject to conditions (34)-( 35) and ( 41) is shown in the lower panel (labeled "VM") of Fig. 2. The interfacial region is evidently close to equilibrium, whereas evaporation occurs in the vdW layer -just like it does under the DIM.The two models are still different, however: the width of the vdW layer in the Vlasov model is comparable to that of the interface, not smaller.This is visible in the lower panel of Fig. 2, but can also be deduced from Eq. ( 40) directly.
To do so, observe that the main contribution to the integral term in (40) comes from the interfacial region (where ρ is large), whereas the contribution of the vdW layer (small ρ) is negligible.Recalling also that the former region is in equilibrium, one can approximate Eq. ( 40) by where ρ (sat) (z) is the equilibrium solution.This equation can be simplified in two different ways.First, recall that ρ (sat) (z) satisfies Eq. ( 42) -hence, the integral term in Eq. ( 47) can be replaced with G(ρ (v.sat) , T ) − Ĝ(ρ (sat) , T ).Second, assume that z is in-side the vdW layer -hence, ρ(z) is small.This allows one to replace µ(ρ, T ) with its low-density limit µ 0 (T ), and Ĝ(ρ, T ) and Ĝ(ρ (sat) , T ), with their low-density asymptotics (16).Eventually, the original integro-differential equation reduces to a differential one, Clearly, the solution ρ(z) has the same spatial scale as ρ (sat) (z) -simply because Eq. ( 48) does not involve other spatial scales.
Equation ( 48) is to be solved with boundary condition (35), rewritten in terms of the relative humidity and saturated vapor density, The boundary condition at the other end is unclear, however -and this is not a technical glitch, but a fundamental issue.It results from the fact that, in the problem under consideration, the neighboring asymptotic zones are on the same spatial scale and, thus, cannot be matched via the van Dyke rule or a similar method.This difficulty can probably be resolved by changing the variables (z, ρ) → (ρ, q = dρ/dz), in which case the interface would correspond to ρ ∼ ρ (l.sat) and the vdW layer, to ρ ∼ ρ (v.sat) .This is, essentially, how matching was handled under the DIM (see Appendix D) -but, in the Vlasov model, the new variables complicate the integral representing the vdW force.
In the end, the following workaround was used.A particular point z m was picked, such that and the vdW layer solutions was 'patched' at z = z m to the interfacial (equilibrium) solution, With ρ (sat) (z) known (pre-computed), boundary-value problem ( 48)-( 49), (51) fully determines the asymptotic solution ρ(z).Note that the 'patching' has been previously used in other problems, and the results have been shown to be asymptotically equivalent to those obtained via matching 24 .Boundary-value problem ( 48)-( 49), ( 51) was solved numerically using the function BVP5C of MATLAB.The patching point z m was chosen such that which automatically satisfies restrictions (50).The computed evaporation rate E can be written in a nondimensional form, relative to in which case E/ ĒV happens to be order one (provided T is not too close to the point where the solution ceases to exist).Typical numerical results are illustrated in the lower panels (labeled "VM") of Fig. 3.The following features should be observed: • The solid curves in panel VM(b) are not extended to small H because the exact boundary-value problem computes much worse than its asymptotic counterpart.
• Panel VM(a) of Fig. 3 illustrates that, for midrange T and H, the asymptotic approach works well (the asymptotic curves in this figure are actually drawn for T ≤ 100 • C, but they are indistinguishable from the exact solution).
• Panel VM(b) illustrates that the asymptotic and exact solutions start to diverge near T = 200 • C.
• Observe that the curves corresponding to different T in panel VM(b) do not collapse onto a single curve [unlike those computed via the DIM and illustrated in panel DIM(b)].In principle, this could be a result of choosing the wrong scale ĒV -and so some other were tested, but none worked.One might conclude that, for the Vlasov model, E depends on T and H in a non-separable way.
• The curves in panel VM(b) of Fig. 3 look like straight lines (but are not ones exactly).
Thus, the DIM and VM both predict an almost linear dependence of E on H, which allows one to compare them to the Hertz-Knudsen Law (HKL).

C. DIM and VM vs. HKL
Rewrite the Hertz-Knudsen Law, given by Eq. ( 2), in terms of the relative humidity To compare this expression to the corresponding results of the DIM and VM, the two latter models should be represented in a similar fashion -say, Once this representation is in place, one can calculate the evaporation/condensation probability, as opposed to just inserting it into the HKL and treating as an adjustable parameter.
One way to obtain a formula for E of form (53) consists in assuming that H → 1 (the vapor is almost saturated) and expanding the DIM or VM solutions in powers of (1 − H).In this expansion, Ē(T ) is the coefficient of the first term.Alternatively, one can determine Ē(T ) by curve fitting, but such an approach is less 'clean' than the asymptotics.
In application to the diffuse-interface model, the expansion in (1 − H) can be found in Ref. 13 , and for the Vlasov model, in Appendix E of the present paper.Under an additional assumption that ρ (v.sat) /ρ (l.sat) ≪ 1, one obtains This expression applies to both DIM and VM, but the resulting Ē depends on which model is used to calculate the profile of the equilibrium interface ρ (sat) (z).
In application to the Vlasov model, formulae ( 53) and (55) are illustrated in Fig. 4. One can see that Ē approximates the slope of the exact curve reasonably accurately.
The dependence of θ on T , calculated via formulae (54)- (55), is illustrated for water in Fig. 5.The following features should be observed: • For T < 100 • C, the DIM and VM both predict that θ is small (in contrast to the classical formulation of the HKL, where θ = 1).
• For T > 250 • C, the DIM and VM both predict that θ is large.This conclusion, however, can be trusted only qualitatively, as ρ (v.sat) /ρ (l.sat) is not necessarily small in this range (making the approximation of viscosity employed for computation of curves (d) and (v) invalid).Note that, even though the coefficient θ in the HKL was initially interpreted as the probability of a molecule to evaporate or condensate, more recent theoretical models (e.g., 6 ) argue that, due to other effects, θ can exceed unity.Thus, the present results are not surprising just because θ > 1, but because θ ≫ 1. • For the most of the temperature range, the DIM noticeably underpredicts θ by comparison with the (more accurate) VM.This is a result of the former's failure in the van der Waals layer.
• It is worth mentioning that the slow-flow approximation [used to reduce the exact momentum equation to Eq. ( 32)] does not impose any restrictions on the results.Indeed, let the Reynolds number be where spatial scale of the interface l F is given by ( 24) and the velocity scale can be expressed through the evaporation rate ĒV of the Vlasov model [given by (52) Estimating expressions ( 56)-( 57) for water, one can show that Re is consistently small for all T between the triple and critical points (reaching the maximum at the latter, where Re ≈ 0.022).

A. Thermodynamics
Following Refs. 15,17, air will be treated as a single fluid with its parameters equal to the 79/21 weighted averages of those of nitrogen and oxygen.Thus, the problem will be formulated for a two-component fluid, representing either a liquid with air dissolved in it, or a mixture of vapor and air.
The thermodynamic state of a multicomponent fluid can be characterized by the temperature T and partial densities ρ i (i = 1 represents the liquid or its vapor, and i = 2 represents the air).Introducing the specific internal energy e(ρ 1 , ρ 2 , T ) and entropy s(ρ 1 , ρ 2 , T ) [satisfying the Gibbs relation (3)], one can define the full pressure and chemical potentials by where is the full density.It can be readily deduced from Eqs. ( 58)-( 60) that This equality is the multicomponent analogue of the pure-fluid identity (6), and Eqs. ( 58)-( 59) are those of Eqs. ( 4)-( 5).The multicomponent Maxwell construction, in turn, is , ρ (a.sat) 2 , T ) = p(ρ G 1 (ρ where a.sat stands for saturated air.To fix the four unknowns -ρ , and ρ (l.sat) 2 -the above three equations should be complimented with the requirement that the saturated air pressure be equal to its atmospheric value, p(ρ Next, the multicomponent version of the Enskog-Vlasov fluid model consists in where c i is the specific heat capacity of the i-th components, R i is its specific gas constant, a ij is the van der Waals coefficient describing the interaction of the i-th and j-th components, and Θ(ρ 1 , ρ 2 ) is the non-ideal part of the entropy.The values of c i and R i can be found in thermodynamics handbooks, and a ij and Θ(ρ 1 , ρ 2 ) should be fitted to the thermodynamic properties of the multicomponent fluid under consideration (See Appendix A 1).
The thermal pressure and thermal chemical potentials are related to the full ones by Using these expressions and identity (61), one can verify that Finally, the low-density asymptotics of the full and thermal chemical potentials are both given by the ideal gas formula, This expression applies only if both ρ 1 and ρ 2 are small, so both fluids can be treated as ideal gases.

B. Governing equations
Let Ψ ij (z) be the one-dimensional potential of the vdW force exerted by component i on component j and vice versa (Ψ ij = Ψ ji ).The van der Waals and Korteweg parameters are now matrices, given by In what follows, the general results will be illustrated using Ψ ij described by formulae ( 26)-( 28) and ( 30), with a and K changed to a ij and K ij (their values for the water-air combination are described in Appendices A 1-A 2).The equations governing slow isothermal dynamics of a binary mixture are where w is the velocity of the mixture as a whole, is the vdW force affecting the i-th component, is the diffusive flux, and the diffusion coefficient D is a known function of ρ 1 , ρ 2 , and T .Observe that J is expressed in terms of the gradients of the chemical potentials, not densities: these two representations are mathematically equivalent, but the former is more convenient in the problem at hand (as well as some others, e.g., Refs. 16,25).
To establish the correspondence between D and the standard diffusivity D which appears in Fick's Law, one needs to adapt the diffusive-flux expression (76) to the ideal-gas limit: set F i = 0 (no vdW force) and replace Ĝi with its small-density asymptotics (70).One should also assume that the density of air exceeds that of vapor (ρ 2 ≫ ρ 1 ) but their gradients are comparable (∂ρ 2 /∂z ∼ ∂ρ 1 /∂z).Eventually, one obtains, to leading order, Comparing this equality to the standard formulation of Fick's Law, one can see that This formula applies only in the low-density limit, which is sufficient for the calculations below.Using identity (69), one can rewrite the momentum equation (74) in the form

C. Boundary conditions
The solution of the governing equations should satisfy calculated together with the full solution.The isobaricity condition (which holds for the multicomponent equations as well) is not sufficient to fix them both.
Boundary conditions ( 79)-( 80) should be complimented with the zero-velocity requirement (36) for the liquid far below the interface, and the zero-viscous-stress requirement (37) for the air far above the interface.

D. Equilibrium interfaces
At equilibrium, there should be no flow (w = 0) and no diffusive flux (J = 0).Substituting expression (75) for F into the momentum equation (78), one can integrate the latter, fix the constant of integration via boundary condition (80) with ρ , and recall Eqs. ( 68) and (71), to obtain Ĝi (ρ were solved numerically using the same algorithm as that for pure fluids.A typical solution is shown in Fig. 6.
The most interesting feature of liquid/air interfaces is the local maximum of the air density.It also arises under the DIM -see Refs. 15,17where it was argued that it emerges because the vdW force exerted by the liquid pulls extra air towards the interface.7][28][29][30] .Since the air density is much smaller than the liquid density, Eqs.(81) can be simplified asymptotically.Firstly, one can neglect the vdW forces exerted by the air on the liquid, its vapor, and itself.Secondly, one can neglect ρ 2 in the expression for G 1 and Ĝ1 (but not in G 2 and Ĝ2 which include ln ρ 2 ).As a result, (81) i=1 reduces to the equation describing the equilibrium liquid/vapor interface in a pure fluid (no air involved), and (81) i=2 becomes Ĝ2 (ρ
For the parameters of water and air at normal conditions, the asymptotic and exact solutions are virtually indistinguishable (as illustrated in Fig. 6).

E. Is diffusion important?
As shown in Ref. 17 , the importance of diffusion is characterized by the following nondimensional parameter: where l is the characteristic interfacial width, ρ is the characteristic density scale, μ is the viscosity scale, etc.
Recall that evaporation of a liquid into its vapor was driven by the van der Waals layer located just outside the interface.The same should be expected for evaporation of liquids into air -hence, one needs to estimate δ specifically for the vdW layer.
To do so, one should use the density and viscosity of air: letting, say, T = 25 • C, one can use Refs. 31,32to obtain (the decimal numbers in the latter formula represent the air's shear and bulk viscosities).For common fluids at normal conditions, the interfacial width is comparable to the scale l F of the vdW force.The latter is given by ( 24) -hence, where with K 11 and a 11 are the Korteweg and vdW parameters of the liquid.For water (see Appendix A), one obtains l ≈ 1 Å, which qualitatively agrees with the spatial scale one can observe in Fig. 6.The diffusion coefficient D, in turn, will be estimated via its low-density asymptotics (77) with ρ 1 = 0.023075 kg m −3 (which is the density of saturated water vapor at 25 • C according to Ref. 19 ) and (which is the diffusivity of water vapor in air at 25 • C according to Ref. 33 ).With these parameter values, expression (83) yields δ ≈ 8.8 × 10 −5 , i.e., the effect of diffusion in the vdW layer is weak.
Further estimates show that diffusion is weak in the interface as well.It is important only at a macroscopic scale, where the diffusive flux matches the evaporative flux emitted by the vdW layer.This part of the setting, however, is not at issue in the present paper.

V. EVAPORATION OF A LIQUID INTO AIR: THE SOLUTION
Let Substituting this ansatz into Eq.( 73), taking into account boundary conditions ( 36) and (79), and omitting the subscript new , one can deduce that Substituting this expression into Eqs.( 78) and (84), one obtains These equations can be viewed as a linear set for the expressions in parentheses on their left-hand sides; solving this set and recalling expressions (75) for F i , one obtains .Curves (v) and (d) are computed using the DIM and VM, respectively.
Since the air density is much smaller than the liquid density, this equation can be simplified the same way the equilibrium problem was: Eq. (85) i=1 can be replaced with its pure-fluid equilibrium version, and in Eq. (85) i=2 , one can set ρ 1 = ρ (sat) 1 and replace Ĝi with its low-density asymptotic (70).
Eq. ( 85) and its asymptotic version were solved numerically, and typical results are illustrated in Fig. 7 together with the corresponding results obtained via the DIM in Ref. 15 .The following features are to be observed: • The dependence of the evaporation rate on the relative humidity is strongly nonlinear (unlike that for evaporation of a pure liquid into its vapor).
• The DIM noticeably underestimates the evaporation rate (similarly to the case of pure liquids).
It is instructive to compare the absolute value of E for the two kinds of evaporation.Using in both cases the Vlasov model, one obtains for water at T = 25 • C and H = 0.5 The huge difference between the two kinds of evaporation is due to the fact that, at 25 • C, air is much denser than water vapor -as a result, the former exerts on the evaporating molecules a much stronger vdW force than the latter.This force pulls the molecules forward -hence, helps evaporation.The back-pulling force exerted by the liquid is the same in both cases, plus it is countered by the pressure gradient.Thus, the net force exerted in the liquid/air system is much more conducive for evaporation than that in the liquid/vapor system.
Note that Eq. ( 85) is much more difficult to solve numerically than its pure-fluid counterpart, for both DIM and VM.The difficulty is probably caused by the presence of two small parameters: the vapor-to-air and vaporto-liquid density ratios -whereas the pure-fluid problem involves only the latter.As a result, it was impossible to extend the curves in Fig. 7 to H ≲ 0.5.For equilibrium interfaces, the difficulties are not as severe (which was the case with pure fluids as well), and there are no difficulties whatsoever for the asymptotic version of the DIM (reduced to an algebraic equation in Ref. 15 ).

VI. HOW CAN THE NEW EFFECT BE OBSERVED/SIMULATED AND WHY HAS THIS NOT HAPPENED ALREADY?
It remains to discuss why the shortcomings of the HKL claimed in the present paper have not been so far noticed by experimental and molecular-dynamics communities.
A. Experiments (i) The author of the present paper found several experimental studies of liquids evaporating into their vapor -but only in a forced setting, where the vapor is sucked out of the container by a pump (e.g., Refs. 34,35).The low pressure created by the pump accelerates the evaporation and makes the predictions of the present paper inapplicable.
More generally, the experimental community do not seem to be concerned with unforced evaporation into vapor, believing that it is similar to the unforced evaporation into air -thus, "why would someone go to significant trouble and expense to do [such] experiments [...] when these can be done in ambient air?" (a private communication from Janet Elliott, Canada Research Chair in Thermodynamics).
(ii) As for liquids evaporating into air, the available measurements of the coefficient θ vary between 0.01 and 1 for the same temperature 4,5 , indicating a problem in the functional dependence where this coefficient appears.
In other words, the experimental community do seem to have noticed the HKL's second shortcoming claimed in this paper.
With this said, experiments with both kinds of evaporation are objectively difficult (hence, potentially inaccurate) because the measurements have to be carried out very near the interface, but without interfering with the evaporative flow.

B. Molecular dynamics
There is a significant body of work where the fluid is approximated by a large set of particles interacting through a potential ϕ(r) involving both repulsive and attractive components.This approach, usually referred to as molecular dynamics, has been applied to evaporation, and recent papers 36,37 claim that the HKL holds with an order-one θ.In the present paper, on the other hand, such is observed only for evaporation of a liquid into its vapor, and only at a mid-range T (see Fig. 5).
Unfortunately, a meaningful comparison between molecular dynamics and Vlasov model is impossible at this stage.
To understand why, note that the choice of the potential ϕ(r) fixes all of the fluid's characteristics, and some of them do not necessarily match the fluid under consideration.The full match can only be accidental, in fact, as none of the commonly-used potentials involves enough adjustable constants to cover the parameter space of a 'real' fluid.As a result, the region where the HKL does not hold could have simply been missed.
Indeed, consider the Lennard-Jones (LJ) potential, used in Ref. 36 with r 0 = 3.41 Å and ϵ = 10.3 meV, to approximate argon.The triple and critical temperatures corresponding to this ϕ(r) are 38 T tr ≈ 79 K, T cr ≈ 158 K, whereas those of the 'real' argon are 39 T tr ≈ 88 K, T cr ≈ 151 K.
Unfortunately, such discrepancies are inevitable, as the LJ potential allows one to explore only a two-dimensional surface (parameterized by ϵ and r 0 ) in the problem's multi dimensional parameter space.Furthermore, Ref. 36 employed a truncated LJ potential (ϕ = 0 for r > 3.2 r 0 ), and the effect of truncation on the fluid's properties is difficult to assess.For example, it can be the reason why liquid/vapor interfaces were observed in Ref. 36 at a lower temperature (T = 76.3K) than both of the above values of T tr .
The mismatch of capillary characteristics is even larger than that of the thermodynamic ones: the LJ value of argon's surface tension -say, at the triple point -is γ tr = 18.6 dyn/cm 38 , whereas its 'real' value is γ tr = 12.6 dyn/cm 40 .
Finally, the vapor-to-liquid density ratio corresponding to the truncated LJ potentials can be very different from that of the 'real' fluid -and this is the most important mismatch of all.
In the simulations of water evaporating into ambient nitrogen reported in Ref. 37 , this parameter was whereas, for 'real' water at, say, 25 • C, the vapor-to-liquid density ratio is smaller by two orders of magnitude, Furthermore, Ref. 37 simulated the case where the vapor and ambient-gas densities were comparable, whereas, in the 'real' atmosphere at 25 • C, this parameter is small, Since the intermolecular force in the vdW layer crucially depends on the density and composition of air, the differences in these characteristics explain the disagreement between the present results and those of Ref. 37 .More generally, to reconcile the Vlasov model and molecular dynamics, one should • either use molecular dynamics with a potential ϕ(r) involving enough adjustable parameters to mimic a common liquid under normal conditions; • or apply the Vlasov model to a fluid whose characteristics match those of the truncated LJ potential.
In the context of the latter approach, note that, for pure water, condition (86) holds if T ≳ 155 • C -and the corresponding values of θ computed via the Vlasov model are order one (see Fig. 5).One can further assume that a small proportion of ambient nitrogen (as in Ref. 37 ) should not alter the VM results too strongly.

VII. SUMMARY AND CONCLUDING REMARKS
This work examines the effect of the van der Waals force on evaporation.The following conclusions have been drawn: (i) For evaporation of a liquid into its vapor, the dependence of the evaporation rate E on the relative humidity H is almost linear -hence, the Hertz-Knudsen Law (HKL) is functionally correct.Yet the evaporation/condensation probability θ, which appears as a coefficient in the HKL, is much smaller than unity, making evaporation much slower than expected.
(ii) For evaporation of a liquid into air, the dependence of E on H is strongly nonlinear, so the HKL does not seem to apply functionally.
In addition to physical conclusions, a technical one has been drawn, which might be important for researchers utilizing the diffuse-interface model: (iii) The DIM fails in a certain region (the vdW layer) at the outskirts of the interface and, as a result, noticeably underestimates the evaporative flux by comparison with the more accurate Vlasov model.
It should be noted, however, that even though the DIM comes short in application to evaporation, it remains to be seen whether it does so in other settings (contact lines, cavitation, etc.).It all depends on whether the DIM solution involves a short-scale boundary layer (vdW layer), making it inapplicable.Furthermore, conclusion (iii) does not apply to a whole class of DIM modelsthose where the density in the interfacial region changes gradually, but the vdW force is not included (e.g., [41][42][43] ).
One should also keep in mind that all conclusions of this work have been drawn using the hydrodynamic approximation of evaporation, which does not describe kinetic effects -such as, for example, the temperature jump associated with the Knudsen layer 3,44 .
To compare the kinetic effects to that of the vdW force, one needs to switch to a kinetic model -e.g., the Enskog-Vlasov equation [45][46][47][48][49][50][51][52][53][54][55] .This is what the author of the present paper initially intended to do, but such a large increase in the model's complexity turned out to be insurmountable in a single stride.
As an alternative to the Enskog-Vlasov kinetic equation, one might use the multi-moment model derived in Ref. 56 .It has a better chance of yielding a relatively simple expression for the evaporative flux, suitable for the use in natural, biological, and industrial applications.

ACKNOWLEDGMENTS
I am grateful to Janet Elliott and Zhi Liang for helping me to understand their work.
Appendix A: The parameters used in the paper This appendix describes how the parameters involved in the DIM and VM can be determined for a specific fluid.The examples considered are water and air; the latter is treated as a mixture of nitrogen and oxygen.

The parameters of the Enskog-Vlasov fluid model
The results described in this subsection were originally reported in Ref. 15,17 and are presented here for completeness.

a. The van der Waals parameter of a pure fluid
To determine the vdW parameter a for a pure fluid, observe that the Enskog-Vlasov expression (9) for the internal energy e(ρ, T ) is linear in ρ.This allows one to determine a as the slope of a linear fit to the empiric dependence of cT − e on ρ, where the heat capacity c is the same as that in the Enskog-Vlasov (kinetic) theoryi.e., 3R for water and 5R/2 for nitrogen and oxygen.For simplicity, the fitting was carried out using only the data on the critical isobar p = p cr , but the resulting straight line fits the isobars p = p cr /2 and p = 2p cr reasonably well too (see Fig. 9(a) of Ref. 17 ).
The values of a determined this way for water, nitrogen, and oxygen are listed in Table II.It also includes the van der Waals parameter of air (calculated as the 79/21 weighted average of those of nitrogen and oxygen, respectively).

b. Thermodynamic properties of a pure fluid
According to the Enskog-Vlasov fluid model ( 10)- (11), the properties of a pure fluid are described by its van der TABLE II.The parameters of H2O, N2, O2, and air: R is the specific gas constant, a is the van der Waals parameter, KD and KV are the values of the Korteweg parameter according to the DIM and VM, respectively.The parameters of air are calculated as the 79/21 weighted averages of the corresponding parameters of nitrogen and oxygen, respectively.
where q (0) ... q (4) are undetermined coefficients and ρ tp is the fluid's density at the triple point (ρ tp is simply a convenient density scale; the fact that, at the triple point, all three phases are in equilibrium is irrelevant).The coefficients q (n) were determined by ensuring that the expressions for p(ρ, T ) and G(ρ, T ) corresponding to (A1) yield the 'correct' -i.e., empiric -values for the critical density, temperature, and pressure, as well as the liquid and vapor densities at the triple point (five equations for the five unknown coefficients).
The values of q (n) for water as determined in Ref. 17 are q (0) = 0.071894, q (1) = 1.4139, q (2) = 8.1126, q ) = −8.3669,q (4) = 4.0238, and those for nitrogen and oxygen can be found in Table 3 of the same paper.The accuracy with which the resulting equations of state approximate the empiric ones is illustrated in Fig. 9(c) of Ref.  c.The van der Waals matrix of a binary mixture When modeling evaporation of water into air, one needs a 11 (water-water interaction), a 22 (air-air interaction), and a 12 (water-air interaction).The first two can be found in Table II, and a 12 can be deduced from a single measurement of the density of air dissolved in water at a certain temperature and pressure.
To do so, consider the equilibrium interface, so that ρ is the density of the air dissolved in water.It generally depends on a 12 -which can, thus, be fixed by fitting ρ (l.sat) 2 to its empiric value.For water at T = 25 • C and p = 1 atm, for example, Ref. 57  Note that this value is specific to the Enskog-Vlasov fluid model used with the Maxwell construction.

The Korteweg parameter
A fluid's thermodynamic properties (discussed above) do not depend on the chosen model of the vdW force, but the capillary properties do.As a result, the DIM and VM correspond to different values of the Korteweg parameter, which will be denoted by K D and K V , respectively.
a.The Korteweg parameter of a pure fluid A fluid's Korteweg parameter can be deduced from a single measurement of its surface tension -say, at the triple point.For water, nitrogen, and oxygen, such measurements can be found in Ref. 58 .
Consider a flat equilibrium interface described by its density profile ρ (sat) (z).Then, according to the DIM, the surface tension is 22 In Ref. 17 , γ D was calculated for water, nitrogen, and oxygen at their respective triple points -and the values of K D were determined, such that (A2) agrees with the corresponding empiric result.These values of K D are listed in Table II.To find K V , one should first derive the Vlasov equivalent of formula (A2).To do so, consider a static macroscopic drop of radius R and calculate the pressure difference between the inside and outside.One should expect it to be of the form γ V /R, where the coefficient γ V is the desired surface tension.
A static (w = 0) density distribution ρ(r) is described by the following reduction of the momentum equation: Write the integral on the right-hand side of (A3) in spherical coordinates (r ′ , β, α) and let the azimuthal angle α be measured from the direction of r.Then, for a spherically symmetric ρ, the integrand does not depend on the polar angle β, and one can reduce (A3) to where (A5) Since the intermolecular potential Φ(r) decays as r → ∞, one can show that the function Ω(r, r ′ ) decays as |r − r ′ | → ∞.One can also verify (see Appendix F) that relationship (81) between Φ and a implies that Impose the following boundary conditions Assuming that ρ(r) describes a spherical drop, introduce the liquid density at its center, and define the drop's radius R by, say, Next, pick R 0 and R ∞ such that R 0 < R < R ∞ and integrate Eq. (A4) from R 0 to R ∞ , which yields where ′2 dr ′ dr, ′2 dr ′ dr, For a macroscopic drop -such that R is much larger than the interfacial width l -these integrals can be simplified.Let R − R 0 and R ∞ − R be much larger than l.Given boundary condition (A7), one can in I 1 set ρ(r) ≈ ρ(r ′ ) ≈ ρ ∞ and obtain Since r ′ ∈ (R ∞ , ∞), the second term in the square brackets is negligible, after which the integral can be estimated via property (A6), In a similar fashion, one obtains Before calculating I 2 , one should symmetrize it with respect to r and r ′ , which yields, after straightforward algebra, The main contribution to this triple integral comes from the region where the solution ρ(r) can be approximated by that for a flat equilibrium interface with its 'midpoint' pinned to r = R, ρ(r) ≈ ρ (sat) (r − R).
Now, take the limit expand I 2 in 1/R, and omit the terms ∼ 1/R 2 and smaller.Changing the variable of integration from α to r ⊥ = √ rr 1 α, one obtains (after straightforward algebra, involving integration by parts) where To ascertain that γ V is the surface tension of the Vlasov model, one should substitute expressions (A10)-(A12) into Eq.(A9) and use (13) to express p through p.The resulting equality shows that the (full) pressure difference between the inside and outside of the drop is indeed γ V /R.Applying expressions (A13)-(A14) [with Φ given by ( 25), ( 27)-( 28), (30)] to water, nitrogen, and oxygen, and making sure that the results match the empiric ones from Ref. 58 , one obtains the values of K V listed in table II.
b.The Korteweg matrix of a binary mixture When modeling evaporation of water into air, one needs K 11 (water-water interaction), K 22 (air-air interaction), and K 12 (water-air interaction).The first two coefficients are listed in Table II, and K 12 is discussed below.
In principle, K 12 could be determined by comparing the characteristics of water/air and water/water-vapor interfaces, but the difference between the two at normal conditions is too small to be reliably measured.Alternatively, K 12 can be deduced from the surface tension of the water/air interface at high pressure, and in Ref. 15 , this was done for the DIM, using the empiric results of Ref. 59 .Unfortunately, the measurements reported in Ref. 59 are scarce (only 6 points) and a bit noisy -thus, the accuracy of the resulting value of K 12 was difficult to assess.
For the present work, the coefficient K 12 was derived from the properties of water/nitrogen and water/oxygen interfaces, both at high pressure, reported in Ref. 60 .The following procedure was used.
Let a flat equilibrium interface in a multicomponent fluid be characterized by ρ (sat) i (z) -say, computed using the Vlasov model.Then, the multicomponent analogue of the Vlasov expression (A13) for the surface tension is where χ ij (z) is given by (A14) with Φ replaced with Φ ij .The coefficient K 12 can be determined by fitting the above theoretical expression to the empiric dependence of γ on p.
An example of such a fit is presented in Fig. 8.One can see that the slope of the theoretical curve is close to that of the empiric one, but the two are separated by a gap.This is because the measurements of Ref. 60 were carried out at 25 • C -whereas, in the present paper, the Vlasov  model is tuned to yield the correct surface tension of pure water at 0 • C. In principle, the gap could be eliminated by retuning the model for 25 • C, but that would only marginally improve the overall accuracy -hence, is not worth implementing.The same procedure was also carried out for the DIM, in which case expression (A15) should be replaced with The calculated values of K 12D and K 12V are listed in Table III.

The viscosity function µ(ρ, T )
When modeling evaporation of water into air, one needs the shear and bulk viscosities, µ s and µ b , of both air and water vapor.The characteristics of liquid water are asymptotically unimportant -as shown in the present work and Ref. 15 for the VM and DIM, respectively (in both cases, provided the liquid's density exceeds those of vapor and air).
In the present work, µ s and µ b of air are calculated using the empiric formulae of Ref. 32 ; and µ s of water vapor, using the IAPWS formulae 61 .
As for µ b of water vapor, there seems to be only one source for it -Ref. 62.In this paper, the results are presented in graphical form, for the interval 58 • C ⪅ T ⪅ 651 • C. The author of the present work digitized them and extrapolated to T = 0 • C. It is worth mentioning here that, at normal conditions, the density of vapor is much smaller than that of air, making the characteristics of the former asymptotically unimportant.
Finally, the effective viscosity µ of the mixture of air and vapor was calculated using the mixture rule proposed in Ref. 63 , Eventually, it was established by trial and error that the proposed algorithm works only if Eq. ( 40) is first integrated with respect to z from −∞ to z ′′ , and then boundary condition (34) The integral on the right-hand side of this equation was evaluated using Simpson's rule [to make the error of the computation consistent with that of the integral on the left-hand side of Eq. (C1)].
Appendix D: Evaporation of a liquid into its vapor: the T → 0 limit of the DIM Evaporation of a pure fluid under the diffuse-interface approximation has been examined in Refs. 14,15using a certain shortcut.In what follows, this shortcut will be reformulated in terms of the standard matched asymptotics, allowing one to estimate the spatial scale of the flow.
The problem will be nodimensionalized using the spatial scale l F of the vdW force given by ( 24), and the following velocity scale w = pl F μ , where p and μ are the pressure and viscosity scales, respectively.Physically, the above choice of w corresponds to the most general regime where the pressure gradient, viscous stress, and vdW force are of the same order 15,17 .Let the density scale ρ be that of liquid and set p = aρ 2 , which reflects the non-ideal part of the pressure [i.e., the second term in the Enskog-Vlasov equation of state (10)].Define the following nondimensional variables: Nondimensionalizing Eq. ( 43) and omitting the subscript nd , one obtains whereas the nondimensional versions of the boundary conditions ( 34)-( 35) look the same as before.One also needs the nondimensional form of the low-density asymptotics ( 12) of p and G, The numerics suggest that the problem involves two asymptotic zones: the interfacial region and van der Waals layer ; the former is near equilibrium and the latter, out of equilibrium.

The vdW layer
In the vdW layer, ρ is small, so one can replace p and G with their low-density expressions (D2).One can also assume µ → 1 as ρ → 0, (D7) which implies that the scale μ used for nondimensionalization is that of the low-density vapor, i.e., μ = µ 0 .Thus, in the vdW layer, Eq. (D1) becomes This equations is to be solved with boundary condition (35).
To simplify Eq. (D8), multiply it by ρ, integrate, use (35) to fix the constant of integration, and change the variables (z, ρ) → (ρ, q), where q is given by (D6).The resulting equation can be written in the form d dρ whereas boundary condition (35) becomes q = 0 at ρ = ρ (v) .(D10) One can readily verify that the solution of Eq. (D9) admits the following asymptotics: where C is an undetermined constant.It can be fixed by matching the inner solution (D11) to asymptotics (D5) of the outer solution.

Matching
The applicability region of the outer (interfacial) solution and that of the inner (vdW layer) solution overlap if ρ (v.sat) ≪ ρ ≪ ρ (l.sat) .
In this region, Eq. (D5) (the inner expansion of the outer solution) should match Eq. (D11) (the outer expansion of the inner solution), which yields C = −1 − ln ρ (v.sat) . (D12) Boundary-value problem (D9)-(D12) determines the function q(ρ) and, more importantly, the dependence of the evaporation rate E on the temperature T and the relative humidity H = ρ (v) /ρ (v.sat) .
T and H can actually be separated by representing E in the form To find the spatial scale of the vdW layer, recall that, nondimensionally, it equals the ratio of the scales of ρ and q in scaling (D14).Recalling also that the spatial variable has been nondimensionalized on l F , one recovers estimate (46).
Eqs. ( 44)-( 45) of the main body of the paper can be recovered by re-dimensionalizing expression (D13).by the equilibrium solution ρ (sat) (z), and the first-order solution exists subject to a certain orthogonality condition which determines E.
The calculation outlined above is straightforward but cumbersome.One can by-pass it by deriving the orthogonality condition directly from the exact boundary-value problem.It should satisfy two requirements: (i) If expanded in ε, the zeroth order of this condition should cancel out.
(ii) The first order should not involve the (unknown) density ρ (l) of the liquid, so that E is the only unknown.
This shortcut was used in Ref. 13 for the DIM -and it can be used, in exactly the same form, for the Vlasov model.
Expanding Eq. (E2) in ε, using the Maxwell construction ( 7)- (8) to ascertain that the zeroth order cancels out, and using identity (6) to simplify the first order, one obtains To simplify Eq. (E3), assume that ρ (l.sat) ≫ ρ (v.sat) -so that the low-density asymptotics (D2) holds for p. Observe also that the largest contribution to the integral on the right-hand side of (E3) comes from the region where ρ (sat) (z) is small -hence, in this region, µ can be replaced with its low-density approximation (D7).
Taking advantage of all these approximations, omitting the O(ε 2 ) terms, and recalling definition (E1) of ε, one can rewrite (E3) in the form where 8,64 where evaporation has been examined using the DIM.The present results suggest that A arises in all hydrodynamic models involving evaporation and vdW force.
Eq. ( 55) of the main body of the paper can be recovered by re-dimensionalizing expression (E5).

FIG. 3 .
FIG.3.The evaporation rate E computed via the diffuse-interface and Vlasov models, scaled by ĒD [see(45)] and ĒV [see(52)], respectively, are presented in the upper and lower panels, respectively.(a) E vs. the temperature T , for two values of the relative humidity H; (b) E vs. H, for two values of T .The solid line shows the numerical solution of the exact equations, the dotted line shows the asymptotic results obtained for the limit ρ (v.sat) /ρ (l.sat) → 0. The temperature in the two left-hand panels ranges from the triple point of water, T ≈ 0 • C, to its critical point, T ≈ 374 • C. The regions where the solution does not exist are shaded (their widths depend on H).

FIG. 4 .
FIG.4.The evaporation rate E computed via the Vlasov model and scaled by ĒV [given by(52)] vs. H, for different values of T .The solid line shows the numerical solution of the exact equations (the same as in panel VM(b) of Fig.3), the dashed line shows the asymptotic results obtained for the limit H → 1.

FIG. 5 .
FIG. 5.The evaporation/condensation probability θ vs. T .Curves (d) and (v) correspond to the DIM and VM, respectively.The straight line marked (h) corresponds to θ = 1, as in the original assumption of Hertz and Knudsen 1,2 .

2 FIG. 6 .
FIG. 6.The equilibrium water/air interface for T = 25 • C. Curves (w) and (a) show the density of water and air, respectively.The exact and asymptotic solutions are shown in solid and dotted line, respectively, but they are virtually indistinguishable.
computed from the pure-fluid problem, one can use this equation to find ρ (sat) 2

FIG. 7 .
FIG. 7. The evaporation rate E vs. H, for T = 25 • C. The solid line shows the numerical solution of the exact equations and the dotted line, the asymptotic result for the limit ρ (a.sat) 2 ≪ ρ (l.sat) 1

FIG. 8 .
FIG.8.The surface tension of water/nitrogen interface vs. the pressure.The solid curve shows the empiric results of Ref.60 , the dashed curve is computed using the Vlasov model.

TABLE I .
Abbreviations used in this paper.
FIG. 1.The pressure p vs. density ρ for the Enskog-Vlasov equation of state

TABLE III .
The nondiagonal Korteweg parameter K12 of water/nitrogen, water/oxygen, and water/air interfaces.(K12) D and (K12) V are calculated according the DIM and VM, respectively.