Investigation of Vortex Structure Modulation by Spume Droplets in the Marine Atmospheric Boundary Layer by Numerical Simulation

Abstract

Direct numerical simulation (DNS) of a droplet-laden, turbulent Couette airflow over a waved water surface is performed modeling the marine atmospheric boundary (MABL) layer carrying idealized spume droplets. Both the instantaneous and mean flow properties, the characteristics of the vortex structures and the momentum exchange between air turbulence and waved water surface and droplet-mediated momentum transfer are investigated. A Eulerian–Lagrangian approach is employed in DNS where full, 3D Navier–Stokes equations for the carrier air are solved in a Eulerian frame, and the trajectories of individual droplets are simultaneously tracked in a Lagrangian frame. The impact of the droplets on the carrier air flow is modeled via a point force approximation. The droplets size is considered in the range of spume droplet sizes observed in MABL. Various water surface roughness and droplet injection scenarios are considered, and both instantaneous and phase-averaged flow fields, the Reynolds stresses and the eigenvalues of the local air velocity gradient tensor are evaluated in DNS. Numerical results show a strong dependence of the droplet-mediated airflow modification on-the-droplet injection mechanism. Droplets injected with the surrounding air velocity effectively mitigate the vortex structures by reducing their swirling strength and suppress the momentum flux from air turbulence to water surface by weakening both ejections and sweeping events, and thus accelerating the mean flow as compared to the droplet-free flow. On the other hand, droplets injected with the velocities of the Lagrangian particles of the water surface enhance both the swirling strength of the vortex structures and air-flow turbulent stresses and decelerate the mean wind. The results also show that these effects of droplet-mediated flow modification become less pronounced as the water surface wave steepness increases.

1. Introduction

Large-scale prognostic models predicting bulk properties of the marine atmospheric boundary layer (MABL) are based on parameterizations of turbulent transfer fluxes of momentum, mass and heat and require detailed knowledge of small-scale exchange processes occurring in the vicinity of the air–sea interface [1]. A turbulent boundary layer developing near the water surface is characterized by the presence of coherent vortex structures, which contribute to the air–water momentum flux and influence the surface drag [2]. Direct measurements of droplet-mediated exchange fluxes in the vicinity of the air–water interface in field conditions and laboratory experiments are technically difficult due to the unsteadiness of the water surface in the presence of wave motions [3]. Investigation of this flow area in the field often becomes even more complicated by air-born sea-spray droplets, which interfere with air–sea interaction and exchange fluxes [4,5,6].

Direct numerical simulation represents an effective research tool and has been successfully employed in the studies of boundary layer flows in the vicinity of waved boundaries, both in droplet-free [7,8,9,10,11] and droplet-laden cases [12,13]. Although the bulk airflow Reynolds numbers in these numerical experiments are typically by orders of magnitude below the values observed in natural MABL, the results provide valuable knowledge of droplet-mediated exchange fluxes, which can be further used in large-scale prognostic models [14]. However, the impact of droplets on the carrier vortex structures has been investigated in detail only for flows over flat (solid) boundaries [15,16]. In these studies, a swirling strength of coherent vortices was explored by evaluating imaginary eigenvalues of the velocity gradient tensor, whereas a quadrant analysis of the velocity field was based on conditional averaging of turbulent velocity fluctuations. These studies show that, under considered conditions and injection scenarios, inertial particles reduce both swirling strength and the number of sweep and ejection events in the boundary layer as compared to the unladen flow case. These findings seem to corroborate an assumption that the impact of sea-spray droplets on MABL may be regarded as a possible reason for the drag reduction observed in hurricanes [6,17]. However, many aspects of the vortex structure modification, related both to the influence of water surface roughness and different possible scenarios of droplets injection, have still remained unexplored.

The objective of this study is to perform direct numerical simulation (DNS) of a turbulent, droplet-laden Couette airflow over a waved water surface. This flow can be regarded as an idealized model of the marine atmospheric boundary layer [7]. Both the instantaneous and mean flow properties and the net impact of the dispersed droplets on the vortex structures and the momentum exchange between air turbulence and water surface are investigated. A Eulerian–Lagrangian approach is adopted where full, 3D Navier–Stokes equations for the carrier air are solved in a Eulerian frame, and the trajectories of individual droplets are simultaneously tracked in a Lagrangian frame. The impact of the droplets on the carrier airflow is modeled via a point force approximation. The droplets diameter is considered in the range of spume droplets sizes observed in the laboratory experiment [18] (typically around 200 micron), which allows us to disregard effects related to droplets deformation. Both instantaneous and phase-averaged flow fields, the Reynolds stresses and the eigenvalues of the local air velocity gradient tensor for two surface wave slope are evaluated for two different droplet injection scenarios. The numerical method is based on an algorithm previously developed for modeling droplet-laden atmospheric boundary layers over progressive surface waves [12] and further adopted below for evaluating the impact of the dispersed spray on MABL vortex structures and air–sea momentum transfer.

2. Governing Equations

The schematic of the problem under consideration is similar to the settings in the previous DNS study [12] (Figure 1). A physical domain with sizes $L_x=6\lambda ,L_y=4\lambda ,L_z=\lambda$, with periodic side boundaries, a solid, flat upper boundary and a waved bottom boundary, $z_s(x,t)$, is considered where the Cartesian coordinates are employed, $0\le x\le L_x;-L_y/2\le y\le L_y/2;z_s\le z\le L_z$. The motion of the air is driven by the upper boundary moving with bulk velocity, $U_0$, and also by the bottom boundary in the form of a prescribed, progressive, two-dimensional (2D) stationary wave of amplitude a, celerity c and length λ propagating in the positive x-direction. The wave celerity is considered sufficiently small (c/U₀ = 0.05, $c/u_{\ast }\approx 1.5$ where $u_{\ast }\approx 0.03U_0$ is the air friction velocity), which is characteristic of “slow” waves as compared to the wind, which is typical under strong wind-forcing conditions [7,8]. In order to avoid coping with a strong geometric nonlinearity caused by the wavy boundary during the integration of the governing equations (to be discussed below), a mapping is introduced transforming the physical domain into a computational domain with a flat bottom boundary, and relating the Cartesian coordinates (x, z) to curvilinear coordinates ($\xi ,\eta$) as:

$$\begin{array}{l}x=\xi -a\exp (-k\eta )\sin k\xi ,\\ z=\eta +a\exp (-k\eta )\cos k\xi ,\end{array}$$

(1)

where $k=2\pi /\lambda$ is the wavenumber. Mapping (1) transforms the wavy boundary at $z=z_s(x,t)$ into a flat-plane boundary at $\eta =0$ (Figure 1). The water surface elevation, $z_s(x,t)$, is defined implicitly by Equation (1) and up to the second order in ka coincides with the Stokes wave solution [19]:

$$z_s(x,t)\approx a\cos k(x-ct)+\frac{a^{2}k}{2}\cos 2k(x-ct)$$

(2)

An additional mapping is also employed for the vertical coordinate in the form:

$$\eta /\lambda =\frac{1}{2}\left(1+\frac{\tanh \tilde{\eta }}{\tanh 1.5}\right),$$

(3)

where $-1.5\le \tilde{\eta }\le 1.5$, so that the grid nodes are clustered in the vicinity of the upper and bottom boundary, at $\eta =0$ ($\tilde{\eta }=-1.5$) and $\eta =\lambda$ ($\tilde{\eta }=1.5$), and stretched towards the middle of the domain (at $\eta =0.5\lambda$, $\tilde{\eta }=0$). The computations are preformed in a reference frame moving with the wave celerity c where the bottom boundary is stationary. Periodic and Dirichlet conditions (to be specified below) are considered at the side and upper and bottom boundary planes, respectively.

A Eulerian–Lagrangian approach is adopted where full, 3D Navier–Stokes equations of the water motion are solved in a Eulerian frame, and the droplets are tracked simultaneously by solving their respective equations of motion in a Lagrangian frame. The Navier–Stokes equations for the carrier air are written in the dimensionless form [1,12,20]:

$$\frac{\partial U_i}{\partial t}+\frac{\partial (U_i{U}_j)}{\partial x_j}=-\frac{\partial P}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}+{\displaystyle \sum _{n=1}^{N_d}{S}_i^n}$$

(4)

where $x_i\equiv (x,y,z)$, $U_j=(U_x,U_y,U_z)$ are the air velocity components, $P$ is the pressure, and $S_i^n$ is a momentum source term (cf. Equation (8) below) contributed by the n-th droplet, and n = 1, …, $N_d$, the latter being the total, constant and number of tracked droplets. Variables in Equation (4) are normalized with velocity scale, U₀, and length scale equal to the wave length, λ. The pressure is normalized with $\rho U_0^2$ where $\rho$ is the air density (≈1.2 kg/m³). The airflow Reynolds number, Re, measures the ratio of the inertial vs. viscous effects and is defined as:

$$\mathrm{Re}=\frac{U_0\lambda }{\nu }$$

(5)

where $\nu$ is the air kinematic viscosity (≈1.5 × 10⁻⁵ m²/s). Equation (4) is supplemented by the incompressibility condition:

$$\frac{\partial U_j}{\partial x_j}=0,$$

(6)

which implicitly defines the pressure field, $P$ (cf. Equation (18) below).

The droplets equations of motion are written as [21]:

$$\frac{d{r}_i^n}{dt}=V_i^n,$$

(7)

$$\frac{d{V}_i^n}{dt}=\frac{f({\mathrm{Re}}_n^{})}{\tau }\left(U_i^n-V_i^n\right)+g{\delta }_{iz}.$$

(8)

In Equations (7) and (8), $r_i^n,V_i^n(i=x,y,z)$ are the n-th droplet coordinate and velocity components, $U_i^n$ is the air velocity at the droplet location, $d/dt$ is the material derivative along the droplet trajectory; g is a dimensionless gravitational acceleration (to be defined below), and ${\delta }_{ij}$ is the Kronecker tensor. The correction in the viscous drag force (accounted for by factor f) is caused by a finite Reynolds number of the droplet and defined as:

$${\mathrm{Re}}_n^{}=\frac{d\left|U^n-V^n\right|}{\nu },$$

(9)

where [21]

$$f\left({\mathrm{Re}}_n\right)=\left(1+0.15{\mathrm{Re}}_n^{0.687}\right).$$

(10)

Dimensionless time scale τ in the drag force on the right hand side of Equation (8) characterizes the ability of the droplet to follow variations of the surrounding air velocity and is defined as:

$$\tau =\frac{d_{}^2}{18\nu }\frac{{\rho }_d}{{\rho }_{}}\frac{U_0^{}}{\lambda }.$$

(11)

Notice that other forces on the droplet by the surrounding carrier air (pressure gradient, added mass, Basset and lift) are neglected as compared to the drag force in the considered case of large droplet-to-air density ratio, ${\rho }_d/\rho \approx {10}^3$, and a sufficiently small droplet diameter (d = 200 µm) [21,22]. Note also that the chosen droplet diameter in DNS is smaller as compared to the viscous (“wall”) length scale (to be discussed below), and the droplet can be treated as a point particle.

The momentum source term in the right hand side of Equation (4), $S_i^n$, is formulated adopting a “point force” approximation [12]:

$$S_i^n=\frac{\pi d_{}^3}{6}\frac{{\rho }_d^{}}{{\rho }_{}}\frac{1}{\tau }\left(V_i^n-U_i(r^n)\right)\left(1+0.15{\mathrm{Re}}_n^{0.687}\right)\frac{w(r^n,r)}{{\mathsf{\Omega }}_g},$$

(12)

where $w(r^n,r)$ is a geometrical weight-factor inversely proportional to the distance between the n-th droplet located at rⁿ = (xⁿ, yⁿ, zⁿ) and the grid node at r = (x, y, z), and Ω_g (r) is the volume of the considered grid cell. Thus, for each individual droplet, eight weight-factors are defined (for each of the surrounding grid nodes) and normalized, so that the sum of partial contributions distributed to these nodes exactly equals the respective total source contribution. Therefore, there is no numerically induced loss or gain of momentum in the droplet–air exchange process (cf., e.g., [12] and references therein for a more detailed discussion).

3. Numerical Method

Since the mapping, Equation (1), is conformal, the following relations between the derivatives hold:

$$\frac{\partial x}{\partial \xi }=\frac{1}{J}\frac{\partial \xi }{\partial x}=\frac{\partial z}{\partial \eta }=\frac{1}{J}\frac{\partial \eta }{\partial z};\frac{\partial x}{\partial \eta }=-\frac{1}{J}\frac{\partial \xi }{\partial z}=\frac{\partial z}{\partial \xi }=-\frac{1}{J}\frac{\partial \eta }{\partial x},$$

(13)

where the Jacobian of the transformation is:

$$J={\left(\frac{\partial \xi }{\partial x_{}}\right)}^2+{\left(\frac{\partial \xi }{\partial z_{}}\right)}^2.$$

(14)

Due to the properties in Equations (13) and (14), the derivatives over the Cartesian coordinates, x and z, can be related to the derivatives over curvilinear coordinates, $\xi ,\eta$, as

$$\begin{array}{l}\frac{\partial }{\partial x_{}}=\frac{\partial {\xi }_{}}{\partial x}\frac{\partial }{\partial {\xi }_{}}+\frac{\partial \eta }{\partial x}\frac{\partial }{\partial \eta }=J\left(\frac{\partial x_{}}{\partial \xi }\frac{\partial }{\partial {\xi }_{}}+\frac{\partial x_{}}{\partial \eta }\frac{\partial }{\partial \eta }\right),\\ \frac{\partial }{\partial z}=\frac{\partial {\xi }_{}}{\partial z}\frac{\partial }{\partial {\xi }_{}}+\frac{\partial \eta }{\partial z}\frac{\partial }{\partial \eta }=J\left(\frac{\partial z}{\partial \xi }\frac{\partial }{\partial {\xi }_{}}+\frac{\partial z}{\partial \eta }\frac{\partial }{\partial \eta }\right).\end{array}$$

(15)

The Laplacian operator is also rewritten as:

$$\frac{{\partial }^2}{\partial x_j\partial x_j}=J\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right)+\frac{{\partial }^2}{\partial y_{}^2}$$

(16)

Equations (4) and (6) for the air velocity are discretized on a staggered grid consisting of $360\times 240\times 180$ nodes, in the $\xi ,y$ and $\tilde{\eta }$ coordinate directions using a second order accuracy finite difference method. The integration of Equation (4) is advanced in time by the second order accuracy Adams–Bashforth method in two stages to calculate the air velocity at each new time step, U_i(t_k+1). First, an intermediate velocity, $U_i^{\ast }$, is computed using the velocity fields at the preceding time steps [23]:

$$U_i^{\ast }=U_i^{}(t_k)+\left(\frac{3}{2}{F}_i^{}(t_k)-\frac{1}{2}{F}_i^{}(t_{k-1})\right)\Delta t,$$

(17)

where the flux, $F_i$, is evaluated as

$$F_i=-\frac{\partial (U_i{U}_j)}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}+{\displaystyle \sum _{n=1}^{N_b}{S}_i^n}.$$

(18)

Further, the new pressure, $P(t_{k+1})$, is computed by solving its respective Poisson equation in the form

$$\frac{{\partial }^{2}P(t_{k+1})}{\partial x_j\partial x_j}=\frac{1}{\Delta t}\frac{\partial U_j^{\ast }}{\partial x_j}.$$

(19)

Equation (19) is solved by iterations by performing; at each iteration step (j) the FFT in the horizontal directions and Gaussian elimination in the vertical direction. The iteration procedure stops when the condition $\left|P_{j+1}-P_j\right|/P_j<0.1\%$ is satisfied. Usually, this condition is met after 3–5 iterations [12]. The new velocity at k + 1 time step satisfying the incompressibility condition (6) is then computed as:

$$U_i^{}(t_{k+1})=U_i^{\ast }-\frac{\partial P{(t_k)}^{}}{\partial x_j}\Delta t.$$

(20)

At the bottom boundary, $\eta =0$, the Dirichlet condition for the air velocity is prescribed to be equal to the surface wave velocity [1] (minus the wave celerity):

$$\begin{array}{l}{U}_x(\xi ,y,0)=c(ka \cos kx(\xi ,0)-1),\\ U_y(\xi ,y,0)=0,\\ U_x(\xi ,y,0)=cka \sin kx(\xi ,0).\end{array}$$

(21)

Periodic conditions for all the fields are prescribed at the side boundaries, and the Dirichlet condition for the air velocity is prescribed at the upper boundary, $\eta =1$, moving with the dimensionless bulk velocity:

$$\begin{array}{l}{U}_x(\xi ,y,1)=1-c,\\ U_y(\xi ,y,1)=0,\\ U_z(\xi ,y,1)=0.\end{array}$$

(22)

The droplet equation of motion is solved by employing the Adams–Bashforth method:

$$V_i^n(t_{k+1})=V_i^n(t_k)+\left(\frac{3}{2}{F}_i^n(t_k)-\frac{1}{2}{F}_i^n(t_{k-1})\right)\Delta t,$$

(23)

$$F_i^n=\frac{f({\mathrm{Re}}_n^{})}{{\tau }_{}}\left(U_i^n-V_i^n\right)-g{\delta }_{iz}.$$

(24)

In Equation (23), the surrounding air velocity at the location of each droplet is obtained via the Hermitian 4th-order accuracy interpolation procedure. The droplet coordinate equation is advanced in time by employing the Adams method as:

$$r_i^n(t_{k+1})=r_i^n(t_k)+0.5\left(V_i^n(t_{k+1})+V_i^n(t_k)\right)\Delta t$$

(25)

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, $\eta$ = 0) or the upper horizontal moving plane (at $\eta =1$) it is re-injected into the flow. Natural and laboratory observations show that spume drops are typically injected in the vicinity of wave crests [18,24,25]. In the present study, the injection model employed in [12] is used: the droplets are injected at heights $0.01<\eta <0.05$ (or $5<\eta \nu /u_{\ast }<25$ in the wall units) above the water surface and randomly distributed at upwind wave slopes in the vicinity of the wave crests. The statistics of spume droplets initial velocity distribution at present are not yet known. Therefore, two different types of injection scenario are considered. Under the first scenario, the droplet velocity at injection equals the surrounding air velocity (the case $V_{inj}=U_{air}$). Under the alternative scenario, the droplet velocity at injection equals the orbital velocity of the water surface particle in an irrotational two-dimensional deep-gravity wave at given x coordinate, Equation (21), ($V_{inj}=U_{surf}$). These two cases of droplet injection are modeled separately.

The airflow bulk Reynolds number in DNS, Equation (5), is set equal to Re = 15,000. The corresponding friction Reynolds number equals ${\mathrm{Re}}_{\ast }=u_{\ast }\lambda /\nu \approx 500$, and is sufficiently large and allows a fully developed turbulent airflow. Experimental observations [25] indicate that for spume droplets with diameters around 200 $\mathsf{\mu }\mathrm{m}$, the ratio of the terminal settling velocity, $V_s^{}=\tau g$ (where $\tau$ is given by Equation (11) vs. the product of the von Kármán constant and the friction velocity, $\kappa u_{\ast }$), is of the order of unity. Thus, the dimensionless gravitational acceleration, g, is prescribed in DNS to satisfy this condition [12].

The air velocity field is initiated as a weakly perturbed, laminar Couette flow. During an initial transient, $0<t<100$, the source terms due to drops on the r. h. s. of Equation (4) is set to zero, and a fully developed turbulent, unladen flow regime sets in. At time $t=100$, droplets are introduced into the flow, with the total (constant) number N_d = 3 × 10⁶, and the initial uniform mass fraction is about 0.038. The droplet volume fraction remains far below 10⁻⁴ and allows neglecting hydrodynamic interactions between neighbouring droplets. The equations of motion of air and droplets, Equations (4)–(8), are solved simultaneously during the time interval $100<t<150$ with the source term “turned off”. During this transient, the droplet dynamics adjusts to the airflow dynamics. At later times, $150<t<200$, the airflow and droplet equations of motion are solved with the source terms “turned on”. During this time interval, a statistically stationary, droplet-laden, two-way coupled flow regime is established.

Similarly, to the previous DNS studies of flows over waved surfaces [7,8,9,10,11,12,13], in the statistical post-processing analysis, phase averaging, equivalent to averaging over an ensemble of turbulent fluctuations, is performed. This averaging (denoted below by angular brackets) is firstly performed over y-coordinate and time t, and further window averaged over ξ—coordinate over six wave lengths as:

$$\langle F\rangle (\xi ,\eta )=\frac{1}{6N_t{N}_y}\sum _{j=1}^{N_y}\sum _{k=1}^{N_t}\sum _{m=0}^{5}F(\xi +m\lambda ,y_j,\lambda \eta ,t_k)$$

(26)

where $F$ is the averaged field, $N_y=240$. Time averaging is preformed over $N_t=500$ realizations in the interval $200<t_k\le 300$. Further, the mean vertical profile of 〈F〉 is obtained by additional averaging along the ξ-coordinate as:

$$\left[F\right](\eta )=\frac{1}{N_x/6}{\displaystyle \sum _{l=1}^{N_x/6}\langle F\rangle ({\xi }_l,\eta )},$$

(27)

where N_x = 360.

The averaging procedure outlined in Equations (26) and (27) is used further to investigate the properties of the air flow vortex structure and sweep and ejections events. In order to define the swirling strength of the vorticity field, a characteristic equation of the velocity gradient tensor, $\Vert \partial U_i/\partial x_j\Vert$, is solved and the complex part of its complex, conjugate eigenvalues is evaluated as in [11,26]:

$${\lambda }_c^2=\frac{3}{4}{\left\{{\left(R+P^{1/2}\right)}^{1/3}-{\left(R-P^{1/2}\right)}^{1/3}\right\}}^2,$$

(28)

$$R=\frac{1}{2}\det \Vert \frac{\partial {\widehat{U}}_i}{\partial x_j}\Vert ,P=\frac{Q^3}{27}+R^2,Q=-\frac{1}{2}\mathrm{tr}{\Vert \frac{\partial {\widehat{U}}_i}{\partial x_j}\Vert }^2,$$

(29)

where

$${\widehat{U}}_i=U_i-\langle U_i\rangle +\left[U_i\right].$$

(30)

It is important to note that the air velocity field (${\widehat{U}}_i$), cf. Equation (30), used for computing ${\lambda }_c^2$ does not include the wave-induced field [11]. Field ${\lambda }_c^2$ is put to zero where P < 0 and averaged according to Equations (26) and (27) to obtain the phase average and mean profiles, $\langle {\lambda }_c^2\rangle$ and $\left[{\lambda }_c^2\right]$.

Sweep and ejections events are characterized by conditional correlations of x- and z- of the velocity fluctuation components (i.e., stresses) defined as [8]:

$$\mathrm{ejection} \mathrm{events}: {\tau }_e=-\langle {U^{\prime }}_x{U^{\prime }}_z\rangle , \mathrm{where} {U^{\prime }}_x^{}=U_x-\langle U_x\rangle < 0; {U^{\prime }}_z^{}=U_z-\langle U_z\rangle > 0$$

(31)

$$\mathrm{sweep} \mathrm{events}: {\tau }_s=-\langle {U^{\prime }}_x{U^{\prime }}_z\rangle , \mathrm{where} {U^{\prime }}_x^{}> 0; {U^{\prime }}_z^{}< 0.$$

(32)

Vertical profiles of the stresses, $[{\tau }_e]$ and $[{\tau }_s]$, are also evaluated according to Equation (27).

4. Results and Discussion

DNS was performed with the same bulk air-flow Reynolds number (Re = 15,000, ${\mathrm{Re}}_{\ast }=u_{\ast }\lambda /\nu \approx 500$, the friction air velocity $u_{\ast }\approx 0.03U_0$), dimensionless wave celerity c = 0.05 and two different wave slopes, ka = 0.1 and 0.2, both for the unladen and droplet-laden cases and different injection scenarios. Monodispersed droplets with diameter d = 200 µm and total number N_d = 3 × 10⁶ were considered. The dimensional time and length scales corresponding to the chosen flow parameters are 0.26 m and 0.3 s; the dimensionless gravity acceleration is g ≈ 0.08 cm/s². The ratio of the dimensionless droplet size over the wall scale equals ${\mathrm{Re}}_{\ast }d/\lambda \approx 0.4$, so the droplets can be treated as point particles.

4.1. Instantaneous Flow Fields

Figure 2, Figure 3 and Figure 4 present instantaneous distributions of the airflow vorticity modulus (in gray scale), $\omega$, defined as:

$$\omega ={\left({\omega }_i{\omega }_i\right)}^{1/2},{\omega }_i={\epsilon }_{ijk}\frac{\partial U_i}{\partial x_j}.$$

(33)

where summation over repeated indices is implied, and droplets respective distributions (symbols) in the central (x,z) and (y,z) planes and in (x,y) plane at the height $\eta \approx 0.04$ above the water surface computed with wave steepness ka = 0.1 for unladen (Figure 2) and droplet-laden cases with two different injection scenarios ($V_{inj}=U_{surf}$, Figure 3, and $V_{inj}=U_{air}$, Figure 4). As illustrated, in the droplet-free case (Figure 2), the instantaneous flow field exhibits numerous boundary layer separation points and occurrences of large-scale vortex structures (also known as “lambda”—vortices) observed also in previous DNS studies of droplet-free boundary layer flows over waved surfaces [7,8,9,27]. The comparison of the vorticity distributions in Figure 2, Figure 3 and Figure 4 shows that droplets injected with the air velocity are confined mostly in the near-surface (z < 0.1) layer and substantially reduce the intensity and quantity of large-scale structures in the droplet-laden air flow as compared to both the unladen-flow and $V_{inj}=U_{surf}$ cases. This impact of droplets on the carrier flow is similar to the effect of attenuation of large-scale vortices by solid particles observed in the vicinity of a flat, solid boundary [14,15]. In the case of droplets injected with the water surface velocity (Figure 4), the air turbulence vorticity attenuation is not so pronounced. The DNS results (not presented here) also show that the vorticity attenuation is also less significant for surface waves with steepness ka = 0.2.

The instantaneous fields in Figure 2, Figure 3 and Figure 4, as well as other characteristics of the flow, are similar to those obtained and discussed in the previous DNS study [12]. So, we refer the interested reader to that discussion, and below we mainly discuss the vortex structure and turbulent momentum modification by the droplets in terms of phase-averaged fields.

4.2. Phase-Averaged Fields

In order to quantify the impact of droplets on the air flow vortex structure and the air–water–surface momentum transfer, the phased-averaged profiles of the complex eigenvalue, ${\lambda }_c^2$, Equation (28), and stresses ${\tau }_e$ and ${\tau }_s$ due to sweep and ejections events, Equations (31) and (32), were evaluated. The results are shown in Figure 5, Figure 6 and Figure 7.

As Figure 5a illustrates, in the unladen-flow case, ${\lambda }_c^2$, which characterizes the swirling strength of vortices, is maximal in the near-surface layer (z < 0.1) above the upwind slope. As the surface roughness increases, the maximum is shifted towards the wave trough (Figure 5b). Droplets injected with water surface velocity ($V_{inj}=U_{surf}$) increase ${\lambda }_c^2$ by about 20% for ka = 0.1 (Figure 5c), and, similarly, reduce ${\lambda }_c^2$ for ka = 0.2, as compared to the droplet-free case. On the other hand, droplets with $V_{inj}=U_{air}$ reduce ${\lambda }_c^2$ to near-zero for ka = 0.1 (Figure 5e) and by about 50% for ka = 0.2 (Figure 5f).

Figure 6 and Figure 7 present distributions of phase-averaged stresses, ${\tau }_e$ and ${\tau }_s$, Equations (31) and (32), obtained for the unladen and droplet-laden flow cases for surface wave steepness ka = 0.1 and 0.2. In the unladen case, the maximum of ${\tau }_e$ is found above the up-wind wave slope (Figure 6a and Figure 7a), whereas the maximum of ${\tau }_s$ is above the windward slope (Figure 6b and Figure 7b). These distributions of ${\tau }_e$ and ${\tau }_s$ indicate that, in the considered case of a “slow” surface wave (as compared to the wind), ejections are dominant at the up-wind slope, and sweeps are mostly found above the windward wave slope.

Distributions in Figure 6c,d show that droplets injected with water surface velocity increase both ${\tau }_e$ and ${\tau }_s$ by about 20% as compared to the unladen flow for ka = 0.1. On the other hand, droplets injected with the surrounding air velocity significantly reduce the stresses (Figure 7e,f). On the other hand, the impact by the droplets in the case ka = 0.2 is less pronounced. This is in accord with the instantaneous flow fields (Figure 3 and Figure 4), which show nearly similar intensity of the vortices in the $V_{inj}=U_{surf}$ case and a drastic reduction in the intensity and the number of the vortices in the airflow in the $V_{inj}=U_{air}$ case.

Figure 8 compares vertical profiles of mean air velocity, $\left[U^{+}\right]$, stresses ($\left[{\tau }_e\right]$ and $\left[{\tau }_s\right]$) and vortex swirling strength, $\left[{\lambda }_c^2\right]$, evaluated from the respective phased-averaged fields as described in Equations (26) and (27), for droplet-free and laden cases with different wave steepness and droplet injection scenarios. The mean air velocity in Figure 8a (ka = 0.1) and Figure 8d (ka = 0.2) is normalized with the friction velocity, $u_{\ast }\approx 0.03U_0$, and the profiles, $\left[U^{+}\right]$, are compared against the viscous asymptotics (in the viscous sublayer region, in a short dashed line) and the Monin–Obukhov (MO) self-similarity law in the logarithmic layer of MABL (in dashed line) [1,17]:

$$U_{visc}={\eta }^{+}\equiv \frac{\eta \nu }{u_{\ast }}, U_{\log }=2.5\ln \left(\frac{\eta }{{\eta }_0}\right),$$

(34)

where ${\eta }_0$ is the roughness length defined as a sum of the viscous length scale and the surface wave amplitude multiplied with a factor defined by a best fit for each considered wave steepness,

$${\eta }_0=0.11\frac{\nu }{u_{\ast }}+\gamma a$$

(35)

where γ = 0.028 for ka = 0.1, and γ = 0.07, for ka = 0.2. Figure 8a,d show that the unladen-flow mean velocity profile is well predicted by the MO theory, Equation (34). In the case ka = 0.1, droplets injected with the surface wave velocity ($V_{inj}=U_{surf}$, blue color) decelerate the air flow (7a,d), and increase both the stresses (8b,e), and the swirl (Figure 8a–c), whereas droplets injected with airflow velocity ($V_{inj}=U_{air}$) accelerate the mean wind, reduce the stresses and the swirling strength.

On the other hand, for wave steepness ka = 0.2 (Figure 8d), the effects of acceleration/deceleration of the wind by the droplets are quantitatively less pronounced, as compared to the case ka = 0.1. The stresses and the swirl are also reduced significantly (by about 30%) by droplets with $V_{inj}=U_{air}$ as compared to the droplet-free flow (Figure 8e,f). However, for this wave steepness, ka = 0.2, the stresses are not enhanced by the droplets injected with the surface wave velocity ($V_{inj}=U_{surf}$) as is observed in simulations with steepness ka = 0.1 (Figure 8b), whereas the swirling strength is reduced by about 30% as compared to the droplet-free flow (Figure 8f).

5. Conclusions

We have performed direct numerical simulation (DNS) of droplet-laden airflow over a waved water surface modeling the marine atmospheric boundary layer (MABL) carrying idealized spume sea-spray droplets. A Eulerian–Lagrangian approach is employed where full, 3D Navier–Stokes equations for the carrier air are solved in a Eulerian frame, and the trajectories of individual droplets are simultaneously tracked in a Lagrangian frame, taking into account the impact of the droplets on the carrier airflow. The droplets diameters are considered close to the typical size of spume droplets in the MABL (around 200 microns). The surface wave shape is prescribed and assumed to be unaffected by either droplets and/or air. Both instantaneous and phase-averaged flow fields, the Reynolds stresses, characterizing the sweeps and ejections events, and the eigenvalues of the local air-velocity gradient tensor are evaluated for two water surface roughness and two different droplet injection scenarios considered in DNS. The specific goal is to investigate the impact of the dispersed spray on MABL vortex structures and air–sea momentum transfer.

The DNS results show a strong dependence of the droplet-mediated modification of the airflow on the spume droplet velocity distribution at injection. Droplets injected with the surrounding air velocity effectively mitigate MABL vortex structures by reducing their swirling strength and suppress the momentum flux from air turbulence to water surface by weakening both ejections and sweeping events, and thus accelerating the mean flow as compared to the droplet-free MABL. On the other hand, droplets injected with the velocities of the Lagrangian particles of the water surface enhance both the swirling strength of MABL vortex structures and air-flow turbulent stresses and decelerate the mean wind. The results also show that these effects of droplet-mediated flow modification become less pronounced as the water surface roughness increases. This duplicity of the simulation results (as far as the droplet-mediated flow modification is concerned) assumes that more yet-unknown details of droplet injection mechanism are to be discovered and incorporated in the numerical code. This however depends on experimental (field and lab) efforts and is thus a subject for future research.