Simulation of Inertial Droplet Dispersion and the Spray Mediated Fluxes in the Atmospheric Boundary Layer Above Waved Water Surface: A Lagrangian Stochastic Model Versus Direct Numerical Simulation

Lagrangian stochastic models (LSM) are widely used to model the dispersion of sea spray droplets injected from the water surface into the marine atmospheric boundary layer (MABL) and for evaluation of the spray impact on the exchange fluxes between the atmosphere and the ocean. While moving through the MABL the droplets pass through the region of high gradients of air velocity, temperature and humidity occurring in the vicinity of the air–water interface. In this case, the applicability of LSMs constructed under the assumption of weakly inhomogeneous flows is questionable. In this work, we develop a Lagrangian stochastic model taking into account the strongly inhomogeneous structure of the airflow in MABL and, in particular, the anisotropy of turbulence dissipation rate. The model constants and the diffusion matrix coefficients are calibrated by comparison of the LSM prediction for the profiles of droplet concentration and the exchange fluxes of sensible and latent heat against the results of direct numerical simulation of turbulent, droplet-laden airflow over a waved water surface.


Introduction
Sea spray is a typical element of the marine atmospheric boundary layer (MABL). At high winds, air-born spume droplets are believed to significantly contribute to the air-sea exchange fluxes of momentum, sensible and latent heat and moisture (see e.g. Veron, 2015 and references therein) the most essential at high winds typical of storm conditions. The most significant contribution to the spray-mediated fluxes is made by large droplets torn by the wind from the crests of waves (see, for example, Kudryavtsev 2006;Kudryavtsev and Makin 2011;Troitskaya et al. 2016Troitskaya et al. , 2018bSroka and Emanuel 2022). One of the possible ways to describe the effect of the droplets on fluxes is based on the use of classical two-phase fluid dynamics equations developed by Barenblatt (1953Barenblatt ( , 1955. In concern with MABL, this approach was applied by, for example, Makin (2005), Kudryavtsev (2006), Bao et al. (2011) and Rastigejev and Suslov (2016), where the idea on the effect of suspended droplets to create a stable stratification of the marine atmospheric boundary layer was exploited. Kudryavtsev and Makin (2011) and Rastigejev and Suslov (2022) proposed a generalization of this approach to take into account the direct effect of momentum exchange between the droplets and the air. In relation to a turbulent boundary layer, this approach requires the use of a number of parameterizations characterized by significant uncertainty. For example, Rastigejev and Suslov (2016) noted that the terms of the heat and moisture balance equations describing the sources associated with droplets contain empirical constants, the significant uncertainty of which is due to the inertia of droplets. Similar problems are noted when describing the spray-mediated source of momentum, which is introduced by Kudryavtsev and Makin (2011). Besides, the equation for concentration based on the description of the balance of deposition and diffusion of particles is well applicable to the description of small particles, but is problematic in the case of large droplets (see, for example, Richter et al. 2019) that make the main contributions to the enthalpy and momentum spray-mediated fluxes.
These problems are completely overcome within the framework of direct numerical simulation (DNS) and large eddy simulation (LES), which, however, have limited practical application. A compromise can be found within the framework of the approach based on the primitive equations of motion and microphysics (Pruppacher and Klett 1987) equations for the droplets in MABL. The main difficulty in this case is associated with the simulation the effect of turbulent fluctuations on the movement and heat exchange of droplets. Lagrangian stochastic models (LSM) are developed to mimic droplets motion in MABL. In the simplest version of LSMs, Markovian chains are used to account for the random velocities, temperatures and masses of droplets and the generalized Langevin equations exploit relaxation models to simulate effects of droplet inertia (see, for example, Edson and Fairall 1994;Edson et al. 1996;Troitskaya et al. 2016). In the LSMs of the second type, the motion of droplets in MABL is described using deterministic equations of motion which include forces that take into account droplets interaction with the surrounding air (e.g. Maxey and Riley 1983) whereas droplet-air heat transfer and droplet evaporation are described in terms of the equations of droplet microphysics (Pruppacher and Klett 1987). A Markovian chain and the stochastic Langevin equations are employed to describe the turbulent fluctuations of the surrounding air velocity, temperature and humidity along the droplet trajectory. This approach was used by Mueller and Veron (2009 to describe sea spray dynamics in MABL. A similar approach was also used for LES of solid particle transport, see for example Glazunov et al. (2016). This approach is also employed in the present work.
The airflow velocity and the scalar fields seen by the droplet moving in a turbulent flow are represented in LSM as a superposition of the ensemble-averaged fields and turbulent fluctuations. For example, the airflow velocity field is represented as: where u is the mean air velocity and u is the velocity fluctuation. There are regular methods to describe the average fields based, for example, on the use of RANS equations with various closure hypotheses, including the case of the atmospheric boundary layer above the waved water surface (e.g., Troitskaya et al. 2014;Makin et al. 1995, etc.). To model the turbulent velocity fluctuations, the Langevin equation expressed in finite differences is usually employed: here a i is the drift term and b i j is the diffusion matrix determined by the properties of the air-flow velocity field. The coefficients a i and b i j in Eq.
(2) are significantly uncertain. There is an approach by Thomson (1987) based on the use of a "well mixing condition" in a flow where turbulent fluctuations are assumed to be Gaussian. In the model by Thompson (1987), b ij is assumed to be a scalar matrix with the equal diagonal elements, namely b ik √ C 0 εδ i j , ε is the rate of dissipation of the kinetic energy of turbulence, C 0 is the Kolmogorov constant, ξ j is the random noise with normal distribution and variance Δt. The classical Kolmogorov constant C 0 6 ± 0.5. (Note however, that the results of LES and LSM of particle dynamics in an atmospheric boundary layer over an urban canopy reported Glazunow et al. (2016) indicate that C 0 is closer to 3-4). There are also models of the diffusion matrix, which reproduce the basic properties of turbulent fluctuations along the trajectories of fluid particles, e.g. the model of the diffusion matrix involving correlated stochastic processes by Moissette et al. (2001), which is constructed to reproduce the second moments of the velocity and scalar fields in the statistically stationary and homogeneous turbulent flow with constant turbulent fluxes.
Given the approximations used by Thomson (1987) in deriving the Langevin equation coefficients it is well applicable in weakly inhomogeneous flows, in which the statistics of velocity fluctuations are close to Gaussian. The approach by Moissette et al. (2001) was also verified by DNS of the turbulent flow with constant second moments of the velocity and scalar fields. These conditions are violated in MABL, since there are large gradients of the average velocity, temperature, humidity and second moments of the turbulent fluctuations in the close proximity to the water surface, where the droplets are injected (see e.g. Monin and Yaglom 2007). Under such conditions, the diffusion matrix is anisotropic and non-diagonal, and the fluctuation statistics differs significantly from Gaussian, as it was shown in numerical experiments (see, for example, Pope 2002;Arcen 2009) in which the Lagrangian statistics of liquid particles in turbulent shear flows were studied. Strictly speaking, in this case, more advanced models of the Lagrangian particle diffusion are required.
At the same time, simpler LSMs used to be applied beyond their applicability. For example, the approach based on a modified model by Thompson (1987) was developed by Veron (2009, 2010) to predict spume droplets dynamics in MABL, however it was not verified. Our goal in this work is to verify the applicability such simplified models to describe the turbulent fluctuations along the particle trajectories and to calibrate the model constants basing on comparison with the results of DNS. The DNS includes solving full 3D Navier-Stoke equations of a turbulent air-flow over a waved water surface and the two-way-coupled Largangian equations of the droplets motion. Although the flow Reynolds number in DNS is by order of magnitude less than the natural analog, it provides valuable information on the details of droplet dynamics which help to evaluate the impact of sea-spray on the momentum, heat and mass exchange occurring in natural MABL (cf. e.g. Peng andRichter 2017, 2019;Druzhinin et al. 2017Druzhinin et al. , 2018. The developed LSM concerns the Langevin equation for the fluctuations in air boundary layer. We considered two versions of the model we have considered two versions of the model, which differ in the drift term and the diffusion matrix. The first version is the modification of the Thompson 1987 equations, which accounts for the anisotropy of the diffusion matrix. The comparison of the model and DNS results also shows that the prediction for the deviation of the trajectories of inertial particles from the trajectories of liquid particles, as well as the isotropic model of turbulent pulsations developed by Mueller and Veron (2009 significantly underestimates turbulent diffusion of the droplets. We also considered an alternative model of the diffusion matrix involving correlated stochastic processes by Moissette et al. (2001), which, in contrast to the Thomson (1987) model, reproduces the second moments of the velocity and scalar fields, including vertical turbulent momentum fluxes. As shown by Moissette et al. (2001), this model reproduces well the dispersion of fluctuations and turbulent fluxes in a turbulent flow with constant velocity shear. Here we also test this model performance in the turbulent boundary layer with strongly inhomegeneous gradients.
The paper is structured as follows. In Sects. 2 and 3, two versions of the LSM and the DNS procedure adapted to the motion and heat and moisture exchange of droplets in MABL are described. In particular, constructing the anisotropic diffusion matrix in the Langevin equation is described. In Sect. 4 the results of DNS and LSM simulations for the average droplet concentrations and turbulent fluxes of momentum, latent and sensible heat to the surrounding air performed for different parameters of the surface waves (slope and celerity) are compared an optimal version and parameters of the model of turbulent fluctuations are proposed.

Lagrangian Stochastic Model of Motion and Heat and Mass Exchange of Heavy Saline Droplets in MABL Over Water Waves
The model described in this section is applied to describe the dispersion of spray droplets in the atmosphere, as well as spray mediated fluxes of momentum and sensible and latent heat fluxes. It extends the model of Troitskaya et al. (2016), to take into account the heat and moisture exchange of droplets with atmospheric air. A new feature of the model is the direct consideration of feedback between droplets and the atmosphere using volume sources of momentum, heat, and moisture associated with spray, which are included in the equations of motion for average thermodynamic and hydrodynamic fields. To describe the scattering of droplets in the atmospheric boundary layer, similarly to Veron (2009, 2010), the Lagrangian equations of individual droplet coordinates r p and velocities V p , and temperatures and masses are used, the coefficients of which, in turn, are determined by average thermodynamic and hydrodynamic fields.
Assuming the droplets small spheres of radius a, temperature T p , and mass m p 4 3 πρ p a 3 , where and neglecting the interaction between the droplets gives the Lagrangian equations are a follows: In Eq. (4), F p 12 Re K (Re)πa 2 ρ a u − V p V p − u is the drag force on the droplet by the surrounding air, u is the instantaneous surrounding airflow velocity. The drop velocity equation, Eq. (4), is written assuming that forces including air-pressure gradient, Basset, added-mass, and shear-induced Saffman's force, can be regarded negligible as compared to the drag and gravity forces for the considered density ratio, ρ p /ρ a ≈ 10 3 (Maxey and Riley 1983). The correction of the drag force K (Re) 1 + 0.15Re 0687 + 0.0175Re 1+4.25·10 4 Re −1.16 in Eq. (4) is introduced as in Clift et al. (1978), taking into account that the droplet Reynlods number, defined as Re 4πκ a T − T p 1 + 0.25Re 0.5 in Eq. (5) is the sensible, diffusive heat flux from the surrounding air to the droplet. F atm→drop is proportional to the product of the instantaneous difference of air temperature at the drop location, T , and drop temperature, T p , and thermal conductivity coefficient, κ'.
H atm→drop 4π D aρ v sat H − H s p 1 + 0.25Re 0.5 in Eq. (6) is the rate of change of the droplet mass. The right hand side of Eq. (6) includes the modified diffusivity of water vapor, D , the saturated vapor density, ρ v sat , and the difference between the surrounding relative humidity and relative humidity at the surface of the drop at saturation, H s p . Expressions for κ and D are taken from Andreas (1989) and cf. also Pruppacher and Klett (1978). These quantities were calculated with the use of a free FORTRAN code (Andreas 2013).
The surrounding hydro-thermodynamic fields of the airflow interacting with the droplet are represented as a sum of three components-mean, wave disturbances and turbulent fluctuations: In the next Subsection, we provide equations for the ensemble-averaged fields and the stochastic equations for the fluctuations.

Equations for Mean Fields
We consider atmospheric boundary layer above the 2D surface wave propagating in xdirection. Then the equations for the turbulence averaged momentum, temperature and humidity are 2D. Similarly to Troitskaya et al. (2016), we formulate the equations in the wave following curvilinear reference frame defined by the mapping: , , t t * , η(x 1 unknown template is a vertical displacement of the water surface due to the wave. In this reference frame the horizontal projection of the RANS equations yields:

Y. Troitskaya
The equations for the turbulence averaged temperature and relative humidity are: and the continuity equation is: In (13), (16) u is the horizontal velocity, w is the vertical velocity,σ 11, 13 are the Reynolds stresses, ρ a is air density, p is pressure, f is a horizontal projection of the density of the bulk force due to the momentum exchange with spray, i.e. momentum transferred from the droplets to the air flow per unit volume per unit time. In (14) and (15), σ T 3 ,σ T 1 , σ q3 and σ q1 are the temperature and humidity fluxes, Q S and Q ρ are the turbulence averaged volume sources of the sensible (Q S ) and latent (Q L ) heat, the amounts of the sensible and latent heat released to the atmosphere from spray droplets per unit time in unit volume. The designation < … > corresponds to statistical averaging over the ensemble of turbulent pulsations. Notice here, that the source terms f , Q S and Q ρ describe the feedback effect of droplets to the average fields in MABL.
In stationary conditions, averaging of (13-15) over x * 1 and integration over x * 3 similarly to Troitskaya et al. 2016 yields the following expression for the conservation of momentum, temperature and humidity fluxes: Integrating Eqs. 18, 19 over x * 3 yields the following expression for the conservation of the momentum, heat and moisture fluxes: where we take into account the boundary conditions are the wave -induced momentum, temperature and moisture fluxes (functions of x * 3 decreasing with the distance from the surface). F x * 3 is momentum delivered from the air flow to spray in the layer from current x * 3 to infinity. If droplets are accelerated by the flow, then F > 0, in the opposite case F < 0. S x * 3 , q x * 3 are sensible heat and moisture released from the spray to the air in the layer above the current x * 3 over the unit area in unit time: In order to obtain the function F(x * 3 ), it is necessary to calculate the horizontal momentum gained by the droplet of the certain radius, when it moves in the air flow above x * 3 and then integrate it over all sizes of droplets ejected from the water surface in unit time over unit square weighted by the size distribution of the ejected spray, i.e. the spray generation function. The functions S x * 3 and q x * 3 are obtained similarly, but the the diffusion heat flux due to the thermal conductivity of the atmosphere and the mass of the water vapor released from the droplet should be calculated.
In order to find the momentum, delivered by the wind to spray in the layer over the distance x * 3 from the water surface, i.e. the function F(x * 3 ) in Eq. (23), it is necessary to calculate the average horizontal momentum gained by the droplets of the initial radius r 0 when they move above x * 3 , p x r 0 , u * , x * 3 solving the system (3-7). Note, that for x * 3 0 it corresponds to the total momentum of the droplets gained during their 'life cycles' from being injected from the water surface to falling down to water. Then, given the number of droplets of a certain radius ejected from the unit water surface in unit time, the spray generation function: In order to find the mass of the water vapor released from the spray to the atmosphere in the layer above the distance x * 3 from the unit water surface in unit time, i.e. the function q x * 3 in Eq. (22), first, it is necessary to calculate the average change of the mass of a single droplet of the initial radius r 0 when it moves above x * 3 , m r 0 , u * , x * 3 solving the system (3-7). Then, given the number of droplets of a certain radius ejected from the unit water surface in unit time, the spray generation function d F dr 0 (r 0 , u * ) yields: The diffusion heat flux from the droplet to the atmosphere is provided by the thermal conductivity and this is opposite to the diffusion flux from the atmosphere to the droplet: F drop→atm −F atm→drop . According to Eq. (5), the sensile heat flux released from a single droplet to the atmosphere is equal to the diffeence of the total heat (or entalpy) and the latent heat (see the details in Troitskaya, et al. 2018b). Then in order to find the sensible heat released from the spray to the atmosphere in the layer above the distance x * 3 from the unit water surface in unit time, i.e. the function S x * 3 in Eq. (21), first, it is necessary to calculate the average change of the total heat (or enthalpy) of a single droplet of the initial radius r 0 when it moves above x * 3 , h r 0 , u * , x * 3 solving the system (3-7). Taking into account, that the average latent heat released to the atmosphere from the single droplet is L v m r 0 , u * , x * 3 and given the number of droplets of a certain radius ejected from the unit water surface in unit time, the spray generation function d F dr 0 (r 0 , u * ) yields finally: Given: where ν t , K T , K q are the eddy viscosity, eddy thermal and moisture conductivity coefficients respectively, the expressions for the average profiles of velocity, temperature and humidity were obtained: The motion and heat and mass transfer of particles in a turbulent boundary layer above waves are greatly affected by turbulent fluctuations in the fields of velocity, temperature, and humidity. Here they are considered within the framework of two versions of LSM. The first one is based on the "well-mixed" condition proposed by Thomson (1987) and in fact repeats the model proposed by Veron (2009, 2010). The second model uses the approach proposed by Moissette et al. (2001), which involves correlated stochastic processes and is designed to reproduce the values of turbulent fluxes and fluctuations in velocity, temperature, and humidity along the trajectories of Lagrangian particles. It should be noted that Veron (2009, 2010) also use the approach of Moissette et al. (2001) to build a model of temperature and humidity fluctuations.

Fluctuation Model Based on the "Well-Mixed" Condition
Turbulent fluctuations of the air flow are described within the Markov chain for velocity, temperature and air humidity fluctuations based on the generalized Langevin equation, respectively: where ξ i are statistically independent Gaussian delta-correlated random processes with variance √ t, the drift term a i −G i j u j is the matrix of inverse relaxation times of the Lagrangian particle velocity turbulent fluctuations, and b i j is the diffusion matrix. To represent the matrix of inverse relaxation times, we used the approach proposed by Thomson (1987), which is based on the "well-mixed" condition in a flow with Gaussian turbulent fluctuations. In this case: hereτ i j u i u j is the stress tensor,τ −1 ik is the inverse stress tensor, U U 0 + u ∼ (x, z, t) is the turbulence averaged velocity vector. Below, this model will be used to describe turbulent fluctuations in a flow in which the mean fields depend only on two components, x and z, and do not depend on time in the reference frame associated with the wave. In this case the stress tensor is as followsτ ⎛ ⎝ τ 11 0 τ 13 0 τ 22 0 And the inverse stress tensor, which enters Eq.
Time derivatives in this case in Eq. (38) turn to zero. Following Thompson 1987 we use here the notation B i j Mueller and Veron 2009). Note, that the classical Kolmogorov's constant C 0 6 ± 0.5, but e.g., Glazunow et al (2016) report C 0 3-4. Besides, as shown by numerical experiments, see for example (Pope 2002;Arcen and Tanière 2009), where the Lagrangian statistics of fluid particles in turbulent shear flows were studied, the diffusion matrix is anisotropic and non-diagonal. This, apparently, reflects the fact that the diffusion of particles in the field of turbulent fluctuations are largely determined by the vortex transport of the energy-carrying scale that have significant anisotropy. Given this we used the expressions for the diffusion matrix, that takes into account the difference in the dissipation rate of the turbulent energy of the components of turbulent flow velocity: where ε i is the dissipation rate of the i-th velocity component. The Kolmogorov constant, C 0 , is the only free parameter of the model. It determines the time scale of the correlation of Markovian random processes describing turbulent fluctuations of the velocity components, temperature and humidity, and, accordingly, the width of the fluctuation power spectra. In this case, the ratio of the correlation scale of the random process and the relaxation time velocity of inertial particles to the equilibrium value determines how much effectively turbulent fluctuations affect the particle velocity. For the matrix of inverse relaxation times of turbulent temperature fluctuations and humidity and diffusion matrices for these values, the same approximations as in Mueller and Veron (2010) were used, based on the approach by Moissette et al. (2001) involving correlated stochastic processes: Mueller and Veron (2010) proposed values for C T 0 C q0 π/2. However, the constants C T 0 , C q0 , along with C 0 , are free parameters, and the sensitivity of the model to them will be studied below.
The solution of the Langevin equations for fluctuations of the velocity, temperature, and humidity components (36), taking into account expressions (38) and (40) for the drift term and the diffusion matrix components, can be written as (see Hall 1975): The model developed by Mueller and Veron (2009, represents the fact that the inertial particles decline from the motions of the elements of the fluid. In the Lagrangian models the separation between the fluid particle trajectory and the heavy particle trajectory is accounted in their autocorrelation coefficients R i (Mueller and Veron, 2009: The components of the separation vector are expressed via t and the components of the velocity vectors of the fluid particle as follows: here However, as will be shown below (see Sect. 4), the modeling of correlation coefficients by expressions (46) leads to a significant underestimation of the effects of particle diffusion in the field of turbulent fluctuations.

Fluctuation Model Involving Correlated Stochastic Processes
Strictly speaking, the model based on the "well-mixed" condition is derived by Thompson (1987) for statistically weakly inhomogeneous and weakly unsteady turbulent flows. In the boundary layer above the wavy surface, the conditions for its applicability are violated. In addition, as mentioned in Mueller and Veron (2010), this model does not reproduce the second moments of the velocity and scalar fields, including vertical turbulent momentum fluxes. The latter, in turn, means a distortion of the phase relation of fluctuations of the velocity components, which noticeably affects the particle trajectory under the action of fluctuations in the flow velocity. Taking into account these circumstances, we have constructed a stochastic model of fluctuations involving correlated stochastic processes, which reproduces the second moments. The principles for constructing such a model are described in the work of Moissette et al. (2001). The generalized Langevin equations for the turbulent fluctuations of the velocity components are as follows: and of temperature and air humidity are as follows: here φ represents fluctuations of temperature, T , or humidity, q . ψ i , ψ φ are five random Gaussian variables with the zero mean value. The variances of ψ i , ψ φ and the pair wise correlators are found from the condition of conservation of all second moments in a statistically homogeneous turbulent flow. The condition of conservation of the fluctuations of the velocity components and turbulent flow dispersion gives the following expression for the correlator: The conservation conditions for the dispersion of temperature and humidity fluctuations, as well as heat and moisture fluxes, give the following relations: In (51) instead of φ either T or q should be substituted. Following the approach proposed by Moissette et al. (2001), we constructed the random variable ψ i as the sum of three statistically independent Gaussian random variables χ k with zero mean and dispersion equal to one: The equations for the coefficients b ik follow from (52). Using the notation introduced earlier u i u j ≡ τ i j they have the following form: System (56) contains 6 different equations in 9 unknowns, i.e. when choosing the coefficients b i k, there is an ambiguity. We propose to use a symmetric matrix in which b i j b ji . Taking into account that in the problem considered here τ 12 τ 23 0, the following system for determining b ik can be proposed: The solution to system (57) can be represented as: We construct the random variable ψ φ in the same way as ψ i , as a superposition of statistically independent Gaussian random variables with dispersion equal to one: The equations for the coefficients follow from (53) and (54) and have the form: The solution of system (59) can be written as: Coefficients b 11 , b 13 , b 33 in (60) are given by the formula (55).
The main goal of this work is to evaluate the applicability of turbulent fluctuation models based on comparison with direct numerical simulation (DNS) data. In this regard, here we will not solve the equations for the mean fields, but take these fields from the DNS, and we will solve only the Langevin Eqs. (35-37) and (50, 51) for the selection of free parameters for which the droplets concentration distribution and mean spray-mediated fluxes of momentum, sensible and latent heat and enthalpy are best reproduced.

Governing Equations and Numerical Method
The schematic of DNS in the present study is similar to that considered by Druzhinin et al. (2018) and Druzhinin (2021) (Fig. 1). A Cartesian framework is employed where the xaxis is oriented along the mean wind, the z-axis is directed vertically upwards, and the y-axis is orthogonal to the mean flow and parallel to the wave front. A two-dimensional, xperiodic surface water wave with amplitude a and wavelength λ is prescribed and considered unaffected by either airflow or droplets. The rectangular computational domain has sizes L x 6λ, L y 4λ and L z λ in the x 1 , x 2 , and x 3 directions, respectively, and the air flow is assumed to be periodic in the x and y directions. DNS is performed in a reference frame moving with the wave phase velocity, c. Thus the bottom boundary representing the wave surface is stationary in this moving reference frame. A no-slip (Dirichlet) boundary condition is prescribed for the air velocity at the bottom boundary, where it coincides with the velocity in the surface wave, and at the top boundary plane moving in the x direction with bulk velocity,U 0 . The air temperature and fractional relative humidity are prescribed as Fig. 1 a Schematic of DNS: L x,y,z are the domain sizes in the horizontal (x), spanwise (y), and vertical (z) directions; a and λ are the surface water wave amplitude and length; T s , f s and T a , f a are the temperature and fractional relative humidity at the water surface and at the top boundary (and their respective profiles in red and blue colour); U 0 is the bulk velocity of the airflow, and g is the acceleration due to gravity T a , f a at the top boundary plane, z L z , and T s , f s , at the water surface. The air and water surface temperatures are set equal to T a 27 • C, T s 28 • C typical of the tropical cyclone conditions. In order to calibrate the constants of the LSM, different surface-wave slopes and surface wave celerities ka 2πa/λ and c/U 0 , are considered in DNS: ka 0.07, c/U 0 0.35; ka 0.1, c/U 0 0.2; ka 0.12, c/U 0 0.15.
In both cases, the relative humidity is prescribed as f a 0.8 and f s 0.98 at the top and surface-wave boundaries, respectively.
The equations for the carrier airflow momentum, temperature and specific humidity are solved in a Eulerian frame simultaneously with the Lagrangian equations of individual droplet coordinates and velocities, and temperatures and masses. The equations for the airflow include the source terms accounting for the contributions individual droplets to the rate of change of air momentum, temperature, and humidity (for details see Druzhinin et al. 2018;Druzhinin 2021). The numerical algorithm in the present study is similar to that employed by Druzhinin et al. (2018) and Druzhinin (2021). The equation are solved on a staggered grid consisting of 360 × 240 × 180 nodes by employing a finite-difference method of second-order accuracy. The integration is performed in curvilinear coordinates (ξ, y, η), which transform the wavy boundary into a plane boundary at η 0. Additional mapping is employed to compress the grid nodes in the vertical direction near boundaries ensuring a sufficient resolution of fine-scale motions. (A more detailed description of the numerical method is provided by Druzhinin et al. 2018 andDruzhinin 2021.) The airflow bulk Reynolds number in DNS is defined as: The bulk air velocity is prescribed as U 0 0.86 m/s and the wavelength is set equal to λ 0.26 m which define the Reynolds number, Re b 15,000. The friction Reynolds number equals Re * u * λ/ν a ≈ 500, where u * is the airflow friction velocity. Note that although Re * is sufficiently large and allows a fully-developed turbulent airflow (as discussed by Druzhinin 2021), the friction velocity in DNS is by an order of magnitude smaller as compared to its natural counterpart. On the other hand, turbulent flow regimes allow us to calibrate the proposed LSM which could be employed further to predict droplets dynamics and its impact on MABL at much higher Reynolds numbers occurring in natural TC conditions. Similarly to the previous studies by Druzhinin et al. (2012Druzhinin et al. ( , 2017Druzhinin et al. ( , 2018 and Druzhinin (2021), the present DNS adopts a reduced-gravity approach in order to adjust droplets the settling velocity to the friction velocity u * corresponding to the considered Re * (which is by an order of magnitude smaller as compared to that in natural conditions). This does not affect the feed-back fluxes but makes the residence times of droplets considered in DNS similar to the adopted estimates (by Andreas 2010) of residence times of spume droplets in air in natural conditions. Accordingly, the dimensional value of the gravitational acceleration was set equal to 8 cm/c 2 .
Statistical sampling of the airflow and droplets is performed during the time interval where the two-way-coupled flow is in a statistically-stationary, turbulent regime (200 < tU 0 /λ < 300). Similarly to the previous studies of airflow over waved surfaces by Sullivan et al. (2000), Yang and Shen (2010), Druzhinin et al. (2012Druzhinin et al. ( , 2017Druzhinin et al. ( , 2018 and Druzhinin (2021), in the statistical post-processing analysis, phase averaging, equivalent to averaging over an ensemble of turbulent fluctuations, is performed over y-coordinate and time, and further window-averaged over ξ -coordinate over six wave lengths, and the mean vertical profiles are obtained by additional averaging along the streamwise coordinate.
The droplet injection mechanism used in the present DNS is similar to the one considered by Druzhinin et al. (2018) and Druzhinin (2021). Accordingly, if a droplet leaves the computational domain via the side boundary plane, it re-enters the domain at the respective opposite side boundary with the same z-coordinate and velocity due to periodicity in x and y directions. If the droplet either reaches the bottom boundary plane (the water surface, η 0) or the upper horizontal moving plane (at η/λ 1) it is re-injected into the flow at random locations distributed at upwind wave slopes in the vicinity of the wave crests with velocities equal (to the first order in ka) to the orbital velocities of the water-surface particles in an irrotational two-dimensional deep-gravity wave, and temperatures equal to the water-surface temperature, T s (cf. Fig. 1).
Since the LSM model presented above has as an input a spray generation function, in DNS we need to evaluate the average number of droplets of a given size injected into the flow per unit time and square of the water surface. Thus the droplets in DNS were distributed over sizes (diameters) into 200 bins between 100 and 300 μm (thus each bin of diameter range d 0 ≤ d ≤ d 0 + d, d 1μm) and the droplets re-injected at each time step in each size bin were counted during the simulation time when the flow was statistically stationary. The resulting size distribution was further normalized with the duration of the considered simulation time interval and the total square of the water surface in DNS. Figure 2 shows the resulting generation function obtained for different DNS runs.
Results of DNS also provided phase-averaged fields of the droplet concentration, C , and droplet-mediated fluxes of momentum and sensible and latent heat, Q V , Q S and Q L , evaluated with the use of the "point force" approximation as follows: where d n is the diameter of n-th droplet, w(r n , r ) is a geometrical weight-factor inversely proportional to the distance between the drop located at r n (x n , y n , z n ) and the grid node located at r (x, y, z), and g (r) is the volume of the considered grid cell (cf. Druzhinin 2021 and references therein for more details). The phased averaged quantities were further averaged over the streamwise coordinate and wavelength to obtain the respective vertical profiles for the fluxes (similarly to the procedure in Sullivan et al. 2000). The respective mean vertical profiles of droplet volume concentration, momentum fluxes, sensible and latent heat, and enthalpy were selected for comparison with LSM.

Comparison of DNS and LSM Results for the Concentration of Droplets in MABL Over Waves
We started by checking how the proposed Lagrangian stochastic models reproduce the second moments of the velocity and temperature fluctuation fields. To do this, we calculated the average values of the dispersion of fluctuations of the velocity, temperature, and air humidity components, as well as the values of the vertical fluxes of momentum, temperature, and humidity on the trajectories of particles binned at the (x 1 , x * 3 ) plane: is the point at the droplet trajectory within a certain bin with the center point x , dx and dz are the grid steps along x 3 , φ-velocity components, temperature or humidity. Figure 3 presents the comparison of the calculation results for two LSMs and DNS. At the same time, calculations for the LSM based on the "well-mixed" condition were carried out with a restriction on the maximum fluctuations of velocity, temperature and humidity fluctuations by 3.5 and 7.5 of the standard deviations of the corresponding values. The Figure  shows that the LSM based on the "well-mixed" condition reproduces well the fluctuations of the horizontal velocity components, as well as temperature and humidity outside the viscous sublayer and somewhat overestimates the vertical velocity fluctuations and underestimates the temperature and humidity fluxes. In this case, the result depends on the given upper limit of fluctuations of thermohydrodynamic values. The quality of reproduction of these LSM values, involving correlated stochastic processes, is comparable. The most significant difference occurs when the models reproduce turbulent shear stress or the vertical momentum flux. The LSM, based on the "well-mixed" condition, fails to reproduce the sign of the vertical momentum flux. This, in turn, means that this version of the LSM incorrectly determines the predominant relationship between the vertical and horizontal velocity phases and, accordingly, the direction of velocity fluctuations. In the presence of a wavy lower boundary, this, in turn, can lead to errors in determining the time the particles stay in the air before colliding with the underlying surface. LSM, involving correlated stochastic processes, is free from this shortcoming. It well reproduces both the magnitude and the sign of the tangential turbulent stress.  Based on a comparison with DNS data, we explored the capabilities of the models for correctly reproducing the particle concentration in the boundary layer and the sources of momentum and heat associated with spray in the equations for mean fields (17,19). For comparison, profiles of integral values averaged over the x * 1 coordinate were used. We calculated the average content of particles in the atmosphere above a certain level G(x * 3 ) (integral of x * 1 -averaged concentration): The profiles of momentum delivered from the air flow to spray,F x * 3 , and sensible heat and moisture, S x * 3 , q x * 3 , released from the spray to the air in the layer above the current z* over the unit area in unit time Eqs. (23-25) were also calculated.
First, the sensitivity of both versions of the LSM to changes in the Kolmogorov constant was studied. Figure 4 shows the profiles G(x * 3 ), F x * 3 S x * 3 and q x * 3 evaluated within the framework of LSM based on the use of the "well-mixed" condition, and LSM, involving correlated stochastic processes, at various values of the Kolmogorov constant C 0 . It is evident that both models are highly sensitive to C 0 when particles are in the turbulent boundary layer. In the viscous and buffer sublayers the profiles G(x * 3 ), F x * 3 S x * 3 and q x * 3 practically coincide.
The   .4 Profiles of the average content of particles in the atmosphere above a certain level G (z*) (a, e), the momentum transferred from the surrounding air to spray,F(z * ) (b, f), sensible heat and moisture, S (z * ) (negative), q (z * ) (positive), (c, g), and the enthalpy (d, h) released from the spray to the air in the layer above the current level z*. The Kolmogorov constant C 0 is equal 2 (black curve), 4 (blue curve), 6 (red curve). a-d are the results of the LSM involving correlated stochastic processes, e-h are the results of the LSM based on "well-mixed" condition. ka 0.07, c/U 0 0.35 of the difference between the profiles of concentration, momentum, sensible and latent heat and enthalpy, evaluated in the framework of DNS and LSM. The norm of the difference was obtained as: Here φ-represents the values of G x * 3 , F x * 3 , S x * 3 , q x * 3 and s x * 3 + q x * 3 , the values z min 0.2, z max 0.8 were chosen in the area of the logarithmic boundary layer. The value of C 0 for ka 0.1, c/U 0 0.2; ka 0.12, c/U 0 0.15 appeared to be 5.01, and for all three cases 5.72. Figure 5 compares the profiles of G x * 3 , F x * 3 , S x * 3 , q x * 3 , and s x * 3 + q x * 3 , evaluated within both versions of the LSM at the optimal value of the Kolmogorov constant  between the prediction of both models and DNS data. However, in our opinion, the use of the LSM version involving correlated stochastic processes is preferable, given that it correctly reproduces the second statistical moments of the turbulent fluctuation fields. It cannot be excluded that a model based on the "well-mixed" condition may work incorrectly in some conditions since it does not reproduce the sign of the turbulent momentum flux in the boundary layer. Note also the model for deviation of the droplet trajectories from the trajectories of carrierair particles, expressed by Eq. (46) sharply underestimates the concentration of particles in the turbulent boundary layer, as well as the droplet-mediated momentum and heat fluxes. In this case, neglecting the effect of deviation of the trajectories gives results that are consistent with the DNS data. This, apparently, reflects the fact that the diffusion of particles in the field of turbulent fluctuations is largely determined by the transfer by energy-carrying eddies. Their characteristic time scale, which can be estimated as the relaxation time T i , is determined by Eq. (47). Figure 6 clearly illustrates that the minimum relaxation time T 3 is an order of magnitude greater than the droplet response time, τ d ρ w ρ a d 2 18ν, even for particles of the maximum size of 300 μm in the present DNS. Therefore, the adaptation of the droplet velocity to the velocity of the energy-containing eddies occurs so quickly that the deviation of the trajectories of the droplets from the trajectories of Lagrangian particles in these eddies can be neglected.
The sensitivity of the calculations results in the framework of both LSM versions to changes in the Kolmogorov constants in equations for scalars was also studied. A change in the constants by a factor of 3 did not lead to insignificant changes in the profiles G (z*), F x * 3 , S x * 3 , q x * 3 and s x * 3 + q x * 3 , which are indistinguishable on the illustrations. Taking this into account, we can recommend the use of the constants proposed in Mueller and Veron (2010).

Discussion and Conclusions
We have presented a Lagrangian stochastic model (LSM) for predicting the dynamics of droplets and the evaluation of droplet-mediated momentum and sensible and latent heat fluxes in the marine atmospheric boundary layer (MABL). The model is based on the equations of motion and heat and moisture exchange of water droplets, which are considered as inertial spherical particles with variable mass and temperature. Note that the proposed model, in addition to the equations of motion and heat and mass transfer of droplets, contains equations for average fields, which include the terms of the spray-mediated fluxes. It shows how we can account for the feedback effect of droplets described within the Lagrangian approach, to the mean fields described within the Eulerian formulation. However, the present work is focused on simulation of turbulent fluctuations at the trajectories of droplets moving in the environment of highly inhomogeneous fields typical for a turbulent boundary layer.
To simulate turbulent fluctuations of the three components of velocity, temperature and humidity, themodel employs a generalized Langevin equation where two types of the matrix of inverse relaxation times are considered. The matrix of the first type was determined from the "well mixing condition" in a flow with Gaussian turbulent fluctuations suggested by Thompson (1987) (also applied by Mueller and Veron 2009, 2014b for MABL). The second type of LSM involving correlated stochastic processes is based on principals suggested by Moissette et al. (2001). The matrix of inverse relaxation times in the Langevin equation is selected to reproduce the variance of fluctuations of three velocity components, temperature and moisture and also the turbulent fluxes, i.e. the second moments of the field of turbulent fluctuations, at the trajectories of the inertial particles. The free parameter of the both models is the Kolmogorov constant.
The relevant model and the Kolmogorov constants have been selected by a comparison with the results of direct numerical simulation (DNS) of a droplet-laden airflow above a waved water surface. The results show that the second model has certain advantages over the first one. First of all, the first model incorrectly reproduces the sign of the shear turbulent stress in the boundary layer. This, in turn, means that the direction of the fluctuation velocity vector in the vertical plane, which provides the diffusion of inertial particles from the underlying surface, is incorrectly reproduced. In addition, it significantly overestimates the dispersion of fluctuations of the velocity components; values close to the DNS data can be obtained only by imposing restrictions on the magnitude of fluctuations. The second model is free from this disadvantage, since the reproduction of the second moments of turbulent fields in it is incorporated in the design of the diffusion matrix. The optimal value of the Kolmogorov constant, C 0 , has been chosen by comparison of the model prediction for the average profiles of content of particles in the atmosphere above a certain level G (z*), the droplet-mediated momentum,F(z * ), and the sensible heat and moisture, S (z * ), q (z * ), released from the spray to the air in the layer above the current z* over the unit area in unit time. This comparison provides C 0 ≈ 5.7, which falls within the range of 4 ± 2 given by Hanna (1981).
The effect of taking into account the deviation of the trajectories of inertial particles from the trajectories of liquid particles was also studied in the framework of the model used in Mueller and Veron (2009, 2014b) (see formula 45). Comparison with the DNS showed a sharp underestimation of the particle concentration in the turbulent boundary layer, and the associated momentum and heat fluxes associated with droplets. At the same time, neglecting this effect gives results that are consistent with the data of a numerical experiment. Taking into account the results obtained, it can be recommended to use LSM involving correlated stochastic processes in determining the spray mediated fluxes, which reproduces the dispersion of field fluctuations and turbulent flows, with the optimal value of the Kolmogorov constant.
Notice, that the LSM is calibrated by the results of a series of DNS at different wave celerity and wave slope. Of course, the question remains open about the applicability of these calibrations at much higher Reynolds numbers (typical for field conditions) than DNS could allow. Although, the turbulent boundary layer regime was realised in the DNS for the flow parameters considered in the work. Due to the use of reduced gravity, it is possible to fulfill the similarity condition with respect to the terminal drop rate and the friction rate between the DNS and high Reynolds number natural conditions. However, the conditions for injection of droplets at high Reynolds number differ significantly from those considered in the DNS. Within the DNS, droplet injection occurs from the buffer region of the turbulent boundary layer. At high Reynolds number corresponding to high wind speeds, the rough flow regime over the water surface is realized (see, for example, Miles 1957). In this case, there is no buffer layer, and the injection of droplets occurs from the region of the turbulent boundary layer. These, and possibly additional factors, can influence the characteristics of turbulent fluctuations along the trajectories of droplets in the air flow, which are modeled by the LSM. The study of such factors goes beyond the scope of this work and requires the use of other computing technologies, for example, the large eddy simulations.