Combined study of evaporation from liquid surface by background oriented schlieren, infrared thermal imaging and numerical simulation

Temperature fields in evaporating liquids are measured by simultaneous use of Background Oriented Schlieren (BOS) technique for the side view and IR thermal imaging for the surface distribution. Good agreement between the two methods is obtained with typical measurement error less than 0.1 K. Two configurations of surface layer are observed: thermocapillary convection state with moving liquid surface and small thermal cells, associated with Marangoni convection, and ”cool skin” with negligible velocity at the surface, larger cells and dramatic increase of velocity within 0.1 mm layer beneath the surface. These configurations are shown to be formed in various liquids (water with various degrees of purification, ethanol, butanol, decane, kerosene, glycerine) depending rather on initial conditions and ambient parameters than on the liquid. Water, which has been considered as the liquid without observable Marangoni convection, actually can exhibit both kinds of behavior during the same experimental run. Evaporation is also studied by means of numerical simulations. Separate problems in air and liquid are considered, with thermal imaging data of surface temperature making the separation possible. It is shown that evaporation rate can be predicted by numerical simulation of the air side with appropriate boundary conditions. Comparison is made with known empirical correlations for Sherwood-Rayleigh relationship. Numerical simulations of water-side problem reveal the issue of velocity boundary conditions at the free surface, determining the structure of surface layer. Flow field similar to observed in the experiments is obtained with special boundary conditions of third kind, presenting a combination of no-slip and surface tension boundary conditions.


Introduction
Evaporation at the free surface of liquid results in cooling of the surface and intensive heat exchange between bulk liquid and atmosphere. A thin surface layer is formed with typical thickness about 1 mm, where temperature drops about 0.5 K. Depending on the conditions, different structures can be observed at the surface and inside the surface layer, including Marangoni convection vortices, cold liquid filaments, motionless "cool skin" and tops of Rayleigh vortices. Despite its small thickness, it is this surface layer which determines the intensity of heat exchange between liquid and atmosphere and evaporation rate. Thus, adequate experimental measurements and theoretical understanding of its thermal structure are crucial for predicting evaporation rates and heat losses in numerous applications involving evaporating liquids. Evaporation from water reservoirs, energy consumption of indoors swimming pool facilities, cooling performance of air-conditioning systems, technological drying processes, geophysical and environmental aspects including life cycle of plankton -just to name a few. Satellite-based remote sensing of the ocean temperature also suffers from ambiguity associated with evaporation: surface temperature, which is actually measured by infrared thermal imaging (IRTI), is different from a e-mail: nickvinn@yandex.ru the bulk temperature, and the difference is determined by evaporation rate.
Average vertical temperature profile below the liquid surface can be approximated by simple exponential fit [1]: where T s is the surface temperature of liquid, T bulk is liquid temperature far from the surface and z is the depth. This profile is governed by two parameters: temperature drop ∆T = T bulk − T s and surface layer thickness δ. Major experimental difficulties arise from the fact that (T bulk − T s ) is usually of order 0.1 K and δ is about 1 mm. This suggests temperature gradients of several hundred K m −1 and thermal fluxes ∼ 10 2 W m −2 . There are two groups of experimental techniques. Methods of the first group allow measuring average values for δ and ∆T . Temperature drop can be found e.g. by simultaneous thermocouple and IRTI measurements [2], and layer thickness is estimated from the total heat flux. Total heat flux is interpreted as molecular one Thus, layer thickness can be found if the total heat flux is known. In laboratory it can be determined from standard thermophysical measurements of liquid sample cooling, or from evaporation rate. In natural conditions small containers are immersed into water reservoir and evaporation rate is measured, or profiler thermoprobes of various types are used [1]. Bulk formulae, relating heat fluxes to temperature and humidity values at various heights, have been elaborated both for laboratory [3] and in situ measurements [4]. These relations allow finding estimates for all heat fluxes constituting the total flux. More information is provided by methods of the second group, which allow measuring 2D fields of temperature and other quantities. These are: shadowgraphy [5], Background Oriented Schlieren (BOS) [11], Particle Image Velocimetry (PIV) performed both in gas and liquid ( [6], [7], [8]), Laser Induced Fluorescence (LIF) [9], and IRTI ( [6], [10]). Principal drawbacks of these methods are well-known. Shadowgraphy requires relatively complex adjustment of experimental setup and the results are mostly qualitative. Nevertheless, it was shadowgraphy to give first hint on complex structure of surface layer [5]. BOS lacks spatial resolution and, like all methods based on refraction, yields distribution, averaged over line-of-sight. PIV also lacks resolution in vertical direction near the surface, and LIF suffers from concentration gradient of emitting particles in the surface layer, which is observationally equivalent to temperature gradient. Also, the techniques with tracer particles immersed into the flow are not truly noncontact. Issues concerning particles behavior near the liquidgas interface are usually limiting the accuracy. IRTI is becoming, with the development of more sensitive devices, one of the major methods of measurement, but it generally observes only very thin layer near the surface since radiation in middle infrared range (3-6 µm) is effectively absorbed by 100 µm of water. Combined usage of PIV and thermal imaging has led to conclusion about complex structure of the flow in upper layer of liquid. It appeared that cold liquid filaments at the surface do not coincide with locations of downward vortical motion [6].
The present study deals with (generally 3D) temperature fields under the surface of evaporating liquids experimentally using BOS for the side view averaged over the tank width and IRTI for top view of the surface temperature, and also by means of numerical simulations. Simultaneous use of two different experimental techniques makes possible the comparison of the results, providing more comprehensive and reliable data and yielding an accuracy estimate. Further comparison with numerical simulations allows making conclusions on the validity of stateof-art models for hydrodynamics of evaporating liquids. IRTI data on surface temperature can provide not only verification of the numerical code, but also boundary condition for separate modeling of air-side and water-side evaporative convection problems. Comparison is also made with known empirical relations for evaporation rate.
Another interesting observation found in literature is the apparent absence of Marangoni convection during water evaporation (see e.g. [12] and references therein). Unlike organic liquids, water does not exhibit small-scale cellular thermal structure flowing along the surface. Instead, a stable film-like layer of motionless liquid with larger thermal cells is observed. This regime, often called "cool skin", will be addressed in the present study too.
The rest of the paper is organized as follows. First, in section 2 the employed experimental techniques and equipment are briefly described. Experimental results are presented in section 3. Numerical models are discussed in section 4 with the results and comparisons given in section 5. Finally, the conclusions are summarized and future perspectives outlined in section 6.

Background Oriented Schlieren
BOS technique was proposed by Meier [11] as a variant of twin-grid schlieren technique complemented with digital cross-correlation processing. Variations of refraction index inside investigated transparent object are measured by digital comparison of two photographs of a background pattern usually merely printed on a sheet of paper. First photograph (reference image) is taken under constant refraction index conditions and the second one (distorted image)through the object being investigated. For evaporation of liquid in tank reference image can be taken with the tank lid on, thereby eliminating evaporation and correspondent heat fluxes. Refraction results in rays deflection and apparent displacement of background pattern elements. The displacement vector field can be found by digital processing similar to PIV. The role of Foucault knife-edge (or a second grid) is thus played by computer. Hence, BOS experimental setup is extremely simple in comparison with conventional schlieren. It merely consists of background pattern, flow under investigation and digital camera. Since the displacement is proportional to refraction index gradient integrated over the camera line-of-sight, it is possible to reconstruct 2D refraction index field (averaged over the tank dimension along line-of-sight) by solving the Poisson equation. Further, density and temperature fields can be derived by solving a system of two algebraic equations: equation of state and Lorentz-Lorenz formula for refraction index -at each interrogation point. Note that this is different from conventional BOS measurements in gases, where Gladstone-Dale relation holds and density can be directly derived from the refraction index. Two other features specific for BOS in liquids are higher (about two orders of magnitude) sensitivity due to larger dn/dT and requirement for flat shape of the tank walls. The latter comes from the fact that the tank with cylindrical (or any curvilinear) walls filled with liquid presents a lens with focal distance determined by the liquid 3D refraction index field. The background pattern image comes out blurred and the deflection angle is dependent on inhomogeneity position inside the tank in this case, so flat walls are mandatory. Also, it should be noted that BOS is sensitive to temperature variations rather than absolute values. In order to obtain absolute values, one has to provide reference temperature value at some point e.g. by measuring bulk liquid temperature with thermocouple. As shown in section 3, this value can be less accurate than BOS results for temperature variations. More information on the optics of BOS measurements and some preliminary results concerning evaporating liquids can be found in [13]. Experiments were conducted in a small glass tank 30 × 50 × 19 mm, in Petri dishes with diameter 9 cm (IRTI measurements only) and in large reservoir 31 × 16 × 25 cm (only for water). Two types of background patterns shown in figure 1 were used: irregular dotted pattern and waveletnoise pattern, proposed in [14]. Dot size in irregular dotted pattern was adjusted for the prescribed distances between the background, water tank and the lens, so that dot image size was about 2-3 pix, which is optimal for cross-correlation interrogation. In contrast, wavelet-noise pattern can be used as universal background for different relative positioning of BOS setup parts, since it contains details of various size. However, it yields slightly larger errors than dotted pattern and is much more vulnerable to image blur. Most of the results were obtained with irregular dotted pattern. Canon EOS 550D SLR camera with 18-55 mm kit lens in 55 mm position was located at 30-40 cm from the tank. Exposure time was 1/20 s for aperture f/14 at ISO 800. Background-to-tank distance was about 1.5 m, resulting in maximal displacements about 5 and 10 pix for water and ethanol, respectively. The camera vertical position was adjusted so that the liquid-gas interface was located in the middle of the frame. This allowed us to minimize the errors associated with multiple ray reflections from the interface. Nevertheless, part of the image below the liquid surface corresponding to about 1 mm of depth had to be cropped out due to heavy blur. Final size of the images was about 1500 × 800 pix. Multi-pass crosscorrelation interrogation with discrete window offset proposed in [15] was employed for displacement evaluation. Typically, three interrogation passes were performed with interrogation window sizes 20×20, 10×10 and 5×5 pix and without window overlapping. IAPWS formulation [16] for water equation of state and empirical dependence of refraction index on density, temperature and wavelength [17] were used to evaluate temperature field in water. Similarly, empirical equation of state [18] and Lorentz-Lorenz formula with volume polarizability 5.41 · 10 −30 m 3 were used for ethanol.
Accuracy of the measurements can be estimated by cross-correlating two reference images. This estimate takes into account the errors of cross-correlation algorithm, lens aberrations, noise of the camera sensor and possible camera displacement. Also, it accounts for refraction index fluctuations which are present even without evaporation. It does not take into account distorted image blur caused by nonlinear refraction index variations [19]. In all cases the estimated error of temperature measurements was about 0.02 K. Bulk liquid temperature was measured by thermocouple with accuracy 0.1 K. Refraction index value correspondent to this temperature and atmospheric pressure was used as boundary condition for Poisson equation at lower boundary. Temperature and relative humidity of the air, which determine the evaporation intensity, were measured by another thermocouple and TESTO-650 gauge.

IR thermal imaging
Measurements were performed with FLIR SC7000-M IR thermal imaging device, operating in wavelength range 2.5-5.5 µm. Since the absorption coefficient of water in this range is at least 80 cm −1 , temperature fields are obtained for surface layer with thickness not more than 100 µm. Image resolution is 640 × 512 pix with sensor pitch 15 µm close to diffraction limit. Camera sensor is cooled down to 80 K with the internal Stirling cooler. The noise equivalent temperature difference is 0.025 K. Two lenses: 25mm f/2.0 and 50mm f/2.0 were used. Camera was tilted down in order to avoid imaging of the camera sensor reflection. However, the incidence angle was small (about 10 • ), hence the surface emissivity was assumed constant, equal 0.96.
Experiments were conducted for different liquids: water, ethanol, butanol, decane, kerosene, glycerine. They can be classified according to their surface tension and volatility (saturated vapor pressure). Glycerine and water possess larger surface tension, whereas ethanol, butanol and decane have similar values. Ethanol is highly volatile, butanol and decane being 9 and 40 times less volatile. Water 01093-p.3 is highly volatile too. These parameters are important for the heat and momentum balance at liquid-gas interface: high volatility intensifies evaporation, thereby increasing the latent heat flux and vertical temperature gradient near the surface. High surface tension, together with temperature variations along the surface, leads to large horizontal velocity differences across the surface layer, leading to "cool skin" regime (see below). Thus, these two parameters, along with ambient and initial conditions, determine the structure of the surface layer formed by evaporation. figure 2 presents typical surface temperature map for ethanol evaporation in Petri dish. Small-scale cells indicating Marangoni convection can be clearly seen. It should be noted that under the thin surface layer thermogravitational Rayleigh convection takes place: the flow field is combination of several large Rayleigh vortices (the exact number being determined the tank aspect ratio) and numerous small Marangoni vortices confined to the surface layer. Liquid at the surface is moving which can be observed either by motion of cellular thermal structures in IR image or by adding powdered coal as tracer particles.

Experimental results and discussion
Surface temperature distribution typical for water evaporation is shown in figure 3. Hot cells are larger. They are surrounded by thin filaments of cold liquid located within the surface layer. In earlier studies these filaments were assumed to be the locus of convective downward motion. However, no correlation was observed between the surface temperature field and liquid velocities beneath the surface [6]. Our experiments confirmed that the surface layer in water exhibits behavior, totally separate from the bulk liquid. Thermal structures at the surface can be completely destroyed by localized IR heating of the surface or by merely covering the tank to cease evaporation. Nevertheless, they are not affected by disturbing under-surface flow with a needle. New filaments can be created using an injector with cold liquid. They are not distorted by fluid motion under the surface. If powdered coal particles are seeded onto the surface, their velocities can be estimated using IRTI video sequence. They are less than 0.1 mm s −1 . However, cold liquid filaments below the surface move with velocity about 1 mm s −1 . Since IR radiation comes from depths not more than 100 µm, a conclusion can be derived that velocity gradient below the surface is very large in comparison with the vorticity of Rayleigh convection vortices. This velocity gradient prevents bulk convective vortices from reaching the surface. Thus, the surface layer presents a thin film of practically motionless liquid ("cool skin"), covering thermogravitational convection vortices. In [20] this behavior was attributed to possible contamination of water surface layer with various surfactants. In fact, small cells typical for Marangoni convection were observed only after special cleaning procedure, applied both to water and tank. Our experiments show that ordinary tap water can exhibit Marangoni convection too. For a given liquid, the regime is determined by the initial vertical temperature gradient. If liquid is hot, evaporation is more intense and Marangoni vortices can break through the motionless liquid film. figure 4 demonstrates that water can exhibit both regimes during the same experimental run. Small Marangoni-like cells are observed in light regions, whereas long cold filaments are found in dark ones. In order to make water exhibit Marangoni convection, it was preheated up to about 60 • C in microwave oven. This increased the vertical temperature gradient and latent heat flux, breaking the motionless cold skin. Note that after the onset of Marangoni convection hot liquid reaches the surface, thereby reducing the temperature gradient. Nonlinear stabilization takes place and small-scale thermal structures disappear. Bidistilled water requires less heating for transition to Marangoni convection: small cells appear even at initial temperature about 40 • C.  Similarly to water, which can be moved to non-typical regime by heating, ethanol can be made to exhibit cool skin. This requires preliminary cooling down to 15 • C. Butanol and decane demonstrate intermediate behavior: at room temperature long filaments are observed, but heating causes transition to Marangoni convection. Glycerine forms cool skin layer up to 120 • C due to small volatility and large surface tension.
The results of simultaneous BOS and IRTI measurements of temperature field for water evaporation in small tank are given in figure 5 (please, refer to electronic version for adequate representation of BOS temperature fields). Comparison of two techniques can be made for temperature profiles near the surface. For that purpose the field obtained by thermal imaging was averaged over the tank width and compared to BOS results at the upper boundary of the domain (intrinsically averaged over BOS line-of-sight). Recall that BOS images were cropped below liquid surface in order to avoid the errors associated with meniscus and multiple reflections of light from the interface. The agreement is very good. BOS temperature is slightly higher, indicating that the upper sublayer about 0.1 mm is not well resolved or is possibly lost during the image crop. However, the difference between two distributions is about 0.05 K, less than accuracy of thermocouple providing the reference temperature value for BOS measurements. Good agreement is related to geometry of the considered flow, which is nearly 2D. If cold filaments are observed near the liquid surface, the complete structure of temperature field is not captured by BOS, since it yields temperature values averaged over line-of-sight.
Analogous comparison for ethanol evaporation is provided in figure 6 for di fferent time intervals after opening the tank cover. As the surface is substantially cooled by evaporation, long cold filament first appears, followed by small cellular structures. Temperature perturbation propagates downwards and can be adequately observed by BOS, even though this technique does not resolve the uppermost

Convection in air
Evaporation problem includes convection on both sides of the interface. Convection in air is affected by two factors: cooling of the near-surface air by water surface, which is cooled by evaporation, and water vapor concentration difference between near-surface and ambient air. These two factors can cooperate or oppose each other e.g. if water surface is cold enough, concentration-induced convection is suppressed by stable thermal stratification. Evaporation rate is thus determined by water surface temperature T w , ambient air temperature T a and relative humidity ϕ 0 far from water surface. Numerous empirical correlations have been proposed, beginning from Carrier's formula in 1918. They can be found in [21] and [22]. Most of them have the form 01093-p.5  where J is evaporation rate, V is the wind velocity, p w and p a are partial vapor pressures at water surface temperature and at ambient air temperature, respectively, and a, b are empirical coefficients. However, this does not include dependence on water reservoir size. Smaller reservoirs exhibit higher evaporation rates, hence correlation Eq. Here D is diffusion coefficient of water vapor, L is the pool size, ρ and η are air density and viscosity, Pr is Prandtl number. Empirical correlations are usually used in numerical simulations in order to provide the mass flux value for the boundary condition at air-water interface (see e.g. [23] and [24]). However, evaporation rate can be predicted numerically, solving Navier-Stokes equations for convection of humid air with specific boundary conditions. First of them comes from the fact that water vapor near the interface is close to saturation. Another one should specify water surface temperature, unless coupled problem is solved for convection both in air and water. It can be either a fixed temperature value (if surface temperature temporal variation is assumed to be slow and one is interested in S h(Ra) relationship) or a local heat balance equation for water surface element where J = −D(∂c/∂y), ∆h is latent heat of vaporization, α is efficient heat transfer coefficient for heat transfer by radiation and through walls of the pool, q l and q r are conductive heat fluxes along the water surface, λ a is thermal conductivity of the air. No-slip boundary condition can be used for air velocity interface since convection velocity in water is typically 10-20 times less than in the air and water level decrease due to evaporation is slow in comparison with convective processes. This results in a closed formulation, which allows one to calculate evaporation rate (and water surface cooling if Eq. (4) is used as boundary condition) without solving the complete coupled problem.

Convection in water
The convective processes in water, forming the observed near-surface thermal layer, are governed by heat and mass exchange at the interface. Evaporation rate at each point is determined by the efficiency, with which the air-side convective flow removes water vapor from the surface. Average evaporation rate (along with initial water temperature) determines the vertical temperature gradient, whereas spatial variations along the surface give birth to Marangoni vortices. Thus, the air-and water-side flows are deeply connected. However, experimental data acquired by IRTI provide an opportunity to separate these problems. Steadystate surface temperature distribution is known and, after appropriate interpolation, can be used as boundary condition for temperature. The crucial issue, which is still left, is boundary condition for velocity. There are several options, which are tested and discussed below in section 5. First, no-slip boundary condition v = 0 can be used. However, it does not account for surface tension and, correspondingly, neglects Marangoni convection. This can be avoided with boundary condition where σ is surface tension, x and y denote horizontal and vertical coordinates, respectively, and the air-side stress is 01093-p.6 neglected. This condition provides large vertical gradient of horizontal velocity, which is observed in experiments. In fact, this boundary condition can be used to estimate vertical gradient of horizontal velocity from experimental data. The right-hand side of Eq. (5) is known from IRTI. Typically, for water it is about 10 s −1 . Expanding the velocity gradient as ∂v x /∂y ∼ ∆v/δ 1 and making use of typical velocity values, observed in PIV or numerical simulations (∆v ∼ v max ∼ 10 −3 m s −1 ), one arrives to an estimate for velocity boundary thickness δ 1 about 100 µm. Bulk convective vortices are characterized by the same velocity difference, but their size is comparable with the tank size, i.e. the vorticity is at least two orders of magnitude less. Note, however, that liquid at the surface is free to move along the surface. It is even forced to do so by Marangoni stress in the right-hand side of Eq. (5). Hence, the regime of motionless cool skin, observed for small vertical temperature gradients can not be obtained using Eq. (5) as boundary condition. It is possible to combine two conditions mentioned above into one condition of third kind where h v is effective momentum transfer coefficient (analogous to heat transfer coefficient) and v 0 is velocity value prescribed at the interface. If there is no wind, v 0 is set to be zero. Relative importance of two terms in the left-hand side of Eq. (6) is adjusted by the value of Bi v = h v L/η -analogue of thermal Biot number. As shown below, for high Bi v values water velocity at the surface is decreased, leading to flow fields resembling those of cool skin regime. An important drawback is that Eq. (6) has no solid physical reasoning and Bi v value is a priori unknown. Another variant is the condition, taking into account the so called surface viscosity η sur f It should be noted that in late 80-s similar problem was considered by Nunez and Sparrow [25] for geometry of vertical open-topped tube. They dealt with a series of numerical models neglecting or taking into account processes of secondary (in comparison with air-side convection) importance: heat exchange between air inside the tube and the ambient, heat conduction through the tube walls, convection in water and radiation heat transfer. However, they did not have experimental IRTI data. Therefore, they did not have means of solving separately water-side problem and overlooked the entire complex of problems associated with surface tension. Also, high Rayleigh numbers could not be addressed at that time because of low computational power.

Computational algorithm and details
Navier-Stokes equations in low Mach approximation are solved. Compressibility effects are taken into account in simplified form, neglecting d p/dt in energy equation. This allows using numerical methods developed for incompressible fluid without considering the sound waves, thus loosening the time step restriction. Additional equation for water vapor transport is solved for air-side problem. Volume condensation is not taken into account, though it can take place locally if water is hot and temperature difference is large. Constant thermal expansion coefficient is assumed for water equation of state in water-side problem or, for large temperature variations, the complete empirical equation of state [16] is incorporated into the model, increasing the computation time. Only 2D simulations are performed. Spalart-Allmaras turbulence model was used, but it appeared that it has minor influence on e.g. evaporation rate because turbulent viscosity near the water surface is close to zero. The geometry of air-side problem presents a lake-like water reservoir surrounded by horizontal flat solid wall, which is assumed adiabatic. Boundary conditions at the water surface are given above, soft boundary conditions are used at upper and lateral boundaries. Ambient air temperature and humidity, which affect evaporation rate, are specified as initial conditions. For water-side problem a rectangular domain with isothermal walls is considered. The upper boundary represents free surface with specified sinusoidal temperature profile T (x) (mimicking IRTI data for cellular surface structures) and boundary conditions for velocity discussed above.
Equations are solved using second-order semi-implicit scheme. Convective terms are approximated with four-points upwind differences. Crank-Nicolson scheme is used for diffusive terms. Time integration is performed by third-order Runge-Kutta method. Pressure is obtained by solving Poisson equation of fractional-step method. For most simulations computational grids 200 × 160 and 100 × 50 are used for air-side and water-side problems, respectively. Grid step is decreasing towards the interface. Code performance was validated by solving the similar problem of natural convection over the heated horizontal plate. The obtained Nusselt-Rayleigh relationship was compared to known empirical correlations [26] Nu = 0.54 Ra 1/4 for laminar flow 0.14 Ra 1/3 for turbulent flow (8) Deviation is small up to Ra ∼ 10 7 , when 3D effects become important. At higher Rayleigh numbers heat exchange is underestimated.

Convection in air
Flow regime is determined by the value of total Rayleigh number, taking into account both the effects of thermal and concentration-induced convection. At negative Rayleigh numbers (cold water case) concentration-induced convection is suppressed by thermal stable stratification. Hydrodynamic velocities are small and humid air is contained within a thin layer near the pool surface. At large positive 01093-p.7 Rayleigh numbers turbulent convection is observed. Evaporation proceeds in chaotic bursts, with plumes of hot humid air emerging sporadically in different locations at water surface. Snapshots and time-averaged fields of relative humidity obtained in simulations for these two regimes are demonstrated in figure 5.1 (please, refer to electronic version). Pool size is 0.614 m, ϕ 0 = 0, T a = 27 • C, T w = 22 • C for case a) and 27 • C for case b). Ra = −5.7 · 10 7 in case a) and 10 8 in case b). White rectangle represents pool surface. Two extreme regimes are separated by the region of laminar convection. Since steady state can be observed for negative and small positive Ra only, calculations were performed until quasi-stationary regime was obtained. Then flow variables and evaporation rate were recorded and averaged over a long time period (e.g., 2000 s for a pool size 1 m). figure 9 presents the comparison between simulation results and known correlations for positive Rayleigh numbers. This implies positive, zero or slightly negative temperature difference between water and air as thermal convection usually dominates over concentration-induced convection. Universal curve is obtained with good accuracy for Sherwood-Rayleigh relationship irrespective of the pool size. Within the range Ra = 10 4 − 10 7 good agreement is observed between simulation results and existing correlations. At higher Rayleigh numbers Sherwood number is underestimated by numerical model and certain discrepancy between different correlations is also found.
Cold water situation when concentration-induced convection is suppressed by thermally stable stratification is less popular. Correlation by Sparrow et al. [27] is the only one known for this case. Nevertheless, this situation is common e.g. in oceanography in upwelling phenomenon. Numerical results are shown in figure 10 as Sherwood number relation to absolute value of Rayleigh number. There is still universal curve in this case, but it lies considerably lower than known correlation. The reason of this discrepancy is unclear and further investigation is needed.

Convection in water
Simulations are performed for a small square tank 2×2 cm. Temperature of the walls is 20 • C. Water surface temperature (except for near-wall boundary layers) varies sinusoidally from 17.5 • C to 18.5 • C with wavelength 0.36 cm, mimicking large cells of warm water and cold water filaments observed experimentally in cool skin regime. Steadystate temperature and horizontal velocity fields are shown in figure 11 for Bi v = 10 and 2 · 10 4 . Distances are normalized with respect to tank size. Note that horizontal velocity is shown only for the uppermost region up to depth 2 mm in order to reveal the structure of the surface layer.
For small values of Bi v number the liquid volume is efficiently cooled with only exception for thermal boundary layers at the walls. Marangoni vortices are more intense than the Rayleigh vortices. Maximal velocities about 3 cm s −1 are observed at water surface. Flow field near the surface, dominated by Marangoni vortices, consists of periodical currents along the surface with downward sinks at locations of temperature minima. The vertical temperature gradient is decreased by warm water reaching the surface. Generally, this scenario corresponds to experimental observations made in Marangoni convection regime (for preheated water or ethanol and other highly volatile liquids). For large Bi v maximal horizontal velocity at the surface decreases down to 0.1 mm s −1 . Marangoni vortices are weak and the flow field of the surface layer is determined by the 01093-p.8 tops of two Rayleigh vortices. The cooled layer is much more thin, except for the central sink of cold water provided by Rayleigh convection. Water at the surface is practically motionless, which corresponds to cool skin regime. No-slip condition and Eq. (5) present extreme cases of boundary condition Eq. (6) for large and small values of Bi v , respectively. No-slip condition does not account for Marangoni vortices, whereas Eq. (5) leads to surface liquid velocities of several cm s −1 , which is incompatible with regime of motionless cool skin. Boundary condition of third kind Eq. (6) yields temperature and velocity fields similar to those observed in experiments, though the value of Bi v number remains an adjustable parameter, which is not derived from properties of liquid and initial (or ambient) conditions. Boundary condition with surface viscosity Eq. (7) is more relevant from physical point of view. Due to presence of velocity second derivative along the surface, it is numerically unstable, which can be fixed by choosing appropriate discretization for this term and performing several iterations to preserve accuracy. The results, similar to shown in figure 11, are obtained for η sur f ≥ 2 · 10 −4 kg s −1 . Note that this value is considerably higher than typically measured in liquid films (10 −5 −10 −7 kg s −1 [28]). The reason of such discrepancy is unclear.

Conclusions
Temperature fields under the surface of evaporating liquids were studied by two independent experimental techniques and also by means of numerical simulations. It was shown that, depending on the liquid surface tension and volatility as well as on initial vertical temperature gradient, two different structures of surface layer can be obtained: motionless cool skin and Marangoni convection. The former is characterized by large cells of warm liquid surrounded by filaments of cold liquid and by surface velocity practically equal to zero. The latter implies small cells and moving liquid surface. Cool skin is typical for liquids with large surface tension and/or small volatility, whereas Marangoni convection usually takes place if surface tension is small or evaporation is intense. Nevertheless, liquids, which typically exhibit cool skin behavior, can be moved into Marangoni mode by simple preheating, providing large vertical temperature gradient and intense evaporation. On the contrary, liquids, which usually form Marangoni small-scale cellular structure, can be cooled down to exhibit cool skin. In particular, water, which, as had been reported, does not exhibit Marangoni convection, can be made to do so. Comparison of the results obtained by BOS and IRTI revealed the lack of spatial resolution of BOS results near the liquid surface. Nevertheless, this technique can be successfully applied to liquid temperature measurements in evaporation problems since the results accuracy is better than 0.1 K. In fact, it is limited by the accuracy of thermocouple, measuring the reference temperature, rather than BOS accuracy. Also, this technique is extremely simple and cheap in realization and, along with IRTI, can be recommended for outdoor measurements.
The performed numerical simulations demonstrate that air-and water-side problems can be separated with help of IRTI surface temperature data. It is shown that evaporation rate can be calculated numerically without making use of empirical calculations. This requires solving air-side problem with appropriate boundary conditions at waterair interface: water vapor saturation and fixed water temperature or water temperature evolution according to local heat balance equation. Evaporation regime is determined by the value of total Rayleigh number, taking into account both the effects of thermal and concentration-induced convection. At high Rayleigh numbers evaporation proceeds chaotically, with plumes of cold humid air formed sporadically at different locations. Universal curves were obtained for Sherwood-Rayleigh relationship in both cases of cold and hot water and compared to known empirical correlations. Simulations of water-side problem and comparison with experimental data revealed the issue of boundary condition for velocity at the surface of evaporating liquid. Boundary condition of third kind was proposed, which enables one to model both variants of the surface layer structure. However, it contains one adjustable parameter, which is not derived from liquid properties and/or initial conditions. Similar situation is observed for boundary condition involving surface viscosity coefficient: simulation provides results being in agreement with experiment only for surface viscosity values considerably higher than found in literature. Future investigations are to solve this problem and to provide physical reasoning for these (or similar) boundary conditions. This might require comprehensive study of near-surface rheology since the surface layer is often observed to exhibit much more viscous behavior than the bulk liquid.
This work was partially supported by Russian Foundation for Basic Research (grant 12-08-01077-a). Authors also acknowledge partial support from M.V. Lomonosov Moscow State University Program of Development.