SEA SPRAY IMPACT ON MOMENTUM EXCHANGES AT THE AIR-SEA INTERFACE DURING HIGH WIND CONDITIONS

At high wind speeds the drag coefficient, characterizing the momentum transfer at the ocean surface, is known to be lower than extrapolation of the existing bulk parameterizations, which were derived at low to medium wind speeds. We hypothesize that sea spray may be responsible for the reduction of the drag, and investigate its effect through direct numerical simulations (DNS). The Lattice Boltzmann method (LBM) is coupled to a Lagrangian particle tracking approach to model numerically the dispersion of sea spray droplets in air turbulence near the air-sea interface during hurricanes and other strong wind events. Our results suggest that the turbulent vortices present near the boundary are damped and/or broken down by the passage of the particles, and that the turbulent Reynolds stress and the production of the turbulent kinetic energy are decreased by the particles. Our results generally agree with previous studies on the turbulence modulation by particles. The streamwise component of the velocity fluctuations is increased, while the spanwise and wall-normal components are decreased compared to the clean channel case. The Reynolds force and the viscous force are reduced and are replaced by the particle feedback force in the force balance. These findings suggest that the sea spray may play an important role in modifying the near surface turbulence during high wind speed events. Our results show that the mean streamwise velocity of the carrier phase is slightly reduced in the logarithmic layer when particles are added to the flow, contrary to the findings in previous studies. Therefore, the particle effect on the drag coefficient remains unclear.


Introduction
In coupled atmosphere-ocean numerical models predicting tropical cyclones (hurricanes), the implementation of accurate boundary conditions at the air-sea interface is critical. The momentum transfer, as well as the heat and moisture fluxes are however difficult to determine in extreme weather conditions, because very few in-situ data exist to corroborate the results of the current numerical simulations. The neutral drag coefficient, which characterizes the momentum transfer between the ocean and the atmosphere, is commonly defined as: where u τ is the friction velocity and is related to the wind stress τ W , τ W = ρ f u 2 τ , ρ f is air density, and U 10 is the neutral wind velocity measured at the 10m reference height and corrected for stability (Foreman and Emeis, 2010;Andreas et al. 2012). Emanuel (1995) theoretically derived a condition in which hurricane intensity would be sustained: the ratio of enthalpy (heat) to momentum surface exchange coefficients C K /C D should be between 1.2 and 1.5. This means that occurrence of intense storms is constrained by how much heat and momentum are exchanged between the ocean and the atmosphere during these events. In particular, this would imply a much lower drag coefficient than the extrapolation of the existing bulk parameterizations, which have been derived based on observations at low-to-moderate wind speeds and are currently used in models. After a series of in-situ observations, Powell et al. (2003) suggested that the drag coefficient is indeed lower than the existing parameterizations at high wind speeds. Since then several hypotheses have been proposed to explain why the drag coefficient is reduced at high winds.
1 Foreman and Emeis (2010) analyzed carefully the expression of the neutral drag coefficient and developed a new parameterization for C D , valid for different wind regimes. At low wind speeds, they observed a nonlinear dependency of the friction velocity u τ on the 10m-height wind speed U 10 . For moderate-to-high wind speeds, they introduced a simpler linear relation: where C m is a modified drag coefficient, and b is a constant with the dimension of a speed. The authors found a good agreement with the existing data from different campaigns. Andreas et al. (2012) pursued the efforts of Foreman and Emeis, and investigated further the relation between the friction velocity and the 10mheight wind speed by extrapolating their reasoning to hurricane-strength winds.
They explained the drag reduction observed in those conditions by considering the wind-wave coupling phenomenon and processing a large amount of data.
Several studies focused on the impact of sea spray on the air-sea exchange coefficients at high wind speeds. Makin (2005) took into account the presence of sea spray during storms to derive a new resistance law (i.e. the expression of the drag coefficient). He did not consider the thermodynamical or the mechanical effect, but instead followed the idea of Barenblatt et al., (2005) about the stratification of the marine boundary layer and the balance of the turbulent kinetic energy. His work was focused on the regime of limiting saturation where the sea-spray droplets dispersed in the air form a dense suspension layer, thus modifying the roughness length and the drag coefficient.
Some studies considered the thermodynamical effect of sea spray on the boundary layer stratification. Andreas and Emanuel (2001) suggested that sea spray droplets could have a cooling effect by giving up sensible heat while being suspended in the air, therefore affecting the intensity of tropical storms. The strat-2 ification of the marine boundary layer was then altered due to the presence of seawater droplets right at the sea surface. More recent works by Kudryavstev et al. (2006Kudryavstev et al. ( , 2011 explored the idea of a thermal impact of the sea spray on the wind turbulence. In the same manner as Andreas et al., the authors suggested that the drag coefficient reduction could be explained through the suppression of the turbulent mixing due to the buoyancy force applied by the spray droplets in the turbulent kinetic energy balance equation. They considered both bubble-formed and spume droplets, the latter being produced from the tearing of breaking wave crests at high wind speeds, and concluded that for spume droplets ejected at the breaking wave height, there could be a major impact on the drag. They parameterized the thermal effect using the Monin-Obukhov similarity theory for stably stratified boundary layers, and noted that their results were in good agreement with Powell's observation regarding the evolution of the drag coefficient for strong wind events. In a similar fashion, Bao et al. (2011) used the Monin-Obukhov framework to investigate the dynamics in the marine boundary layer. The authors described how sea spray droplets exchange sensible heat with the atmosphere and enhance buoyancy at the sea surface. They also explored the mechanical effects of sea spray, and showed that the mechanical impact counterbalances the thermal effect by reducing the friction velocity, and stabilizing the marine boundary layer. They explained how the presence of sea spray intensifies the storms by decreasing the drag through numerical simulations.
Most recently,  performed a series of direct numerical simulations to highlight the role of sea spray in the momentum transfer between the atmosphere and the ocean. The sea spray droplets were modeled as heavy inertial particles, suspended in a turbulent airflow. They used a common 3 approach from CFD (Computational Fluid Dynamics) and Engineering, the Eulerian-Lagrangian approach, where the droplets were represented by solid pointwise particles whose individual trajectory was tracked over time in a Lagrangian way. At the same time, the fluid governing equations (the traditional Navier-Stokes equations plus a feedback term from the sea spray) were solved in an Eulerian framework. They discussed the fact that for typical diameters of sea spray droplets (typically from 10m to 1mm), the mechanical effect would dominate over the thermal effect. By studying the momentum budget, Richter and Sullivan introduced a "spray" stress, which compensates the decrease in Reynolds stress by providing a feedback effect to the turbulence. They explained that the drag coefficient based on the total stress remains almost unchanged in the presence of sea spray, while the drag coefficient based on the turbulent stress is reduced.
Within the framework of this project, we follow the approach of . Our work focuses on the impact of an idealized sea spray on the atmospheric turbulence in the marine boundary layer, using the Lattice Boltzmann Method (LBM) to determine what effects sea spray droplets have on the evolution of the drag coefficient in high wind conditions. We model this problem by simulating a turbulent channel flow where particles are introduced in place of the sea spray droplets to reproduce their effects on the flow.
The approach taken here represents an example of turbulent particle-laden flows, which constitute an entire branch of CFD. More specifically, turbulence modulation by particles (droplets, bubbles or solid particles) forms the main interest of this branch. The study of particle-laden turbulent flows has been ongoing for the last twenty-five years. However, because the full resolution of such flows is computationally very demanding, a hybrid approach has been commonly 4 adopted by many scientists, where the particles are described in a Lagrangian way while the fluid is solved in an Eulerian framework. This is called the Lagrangian point-particle approach (Balachandar and Eaton, 2010), and it enables to track the particle positions by solving the equations of motion in the Lagrangian framework, providing that the particles are small compared to the characteristic length scale of the turbulence. In addition, the mass, the momentum and the energy of the dispersed phase are also being solved for the particles in the Lagrangian framework. In a similar fashion as Richter and Sullivan (2013), the particle-tracking method is retained in our study to simulate sea spray droplets without fully solving the flow around each individual droplet, hence saving great amount of computational effort and time.
The behaviour of dispersed multiphase flows can be classified according to the interactions between the carrier phase and the dispersed phase. Elghobashi (Elghobashi, 1994) proposed a regime map based on the interactions between the fluid phase and the dispersed phase. The original map, drawn for homogeneous turbulence by Elghobashi is presented here.
The figure 1 illustrates the different interactions that can be observed in such flows in terms of the two parameters, volume fraction of particles Φ p , and the particle response time τ p normalized by the turbulence time scale τ e or τ K . Here, τ K is the Kolmogorov time scale, given by (ν/ε), where ν is the kinematic viscosity of the fluid and ε is the dissipation rate of the turbulent kinetic energy. When Φ p is below 10 −6 , the flow affects the particle trajectories but the dispersed phase has a negligible effect on the turbulence: the interaction between the two phases is called one-way coupling (i.e. from fluid to particles). For intermediate concen- trations ranging from 10 −6 to 10 −3 , the momentum transfer from the particles to Figure 1: Original regime map existing between the carrier phase and the dispersed phase (Elghobashi, 1994) the turbulence is large enough to alter the turbulence structure (Elghobashi, 1994) and the interaction is described as two-way coupling i.e. (from fluid to particles, and from particles to fluid). In the case of dense suspensions, when Φ p is above 10 −3 , interactions between particles become possible, such as collisions, and the interaction is called four-way coupling (i.e. between carrier fluid and particles, and among particles).
On one hand, one-way and two-way coupled systems have been extensively investigated, both numerically and experimentally (Balachandar and Eaton, 2010). On the other hand, the four-way coupling implies that the dynamics of droplets should be taken into account, such as breaking up or coalescence, and 6 represents a computational challenge. Hence, few studies have been devoted to this topic. Although this approach would be more realistic in the sea spray problem, it also requires a much more complex numerical treatment. We therefore concentrate our efforts on modeling the two-way coupling between sea spray droplets and the air flow, assuming the latter can be represented as pointwise particles dispersed in a turbulent bounded flow. Hence, the interactions among particles are neglected in this study.
This dissertation is structured as follows: the introduction of the governing equations for the dispersed phase and the fluid, as well as the main assumptions on which our model is based, will be the focus of chapter 2. The Eulerian-Lagrangian method and the Lattice Boltzmann method will be described in chapter 3, and the results of our several DNS runs will be presented in chapter 4. We will conclude with discussions of the main results and the general conclusions of our study.

Mathematical Description
In this study the impact of sea spray on the air turbulence in high wind conditions is investigated. A series of DNS is performed to model the interactions between the seawater droplets and the turbulent airflow at the base of the marine boundary layer. The problem is schematized in a rather simple manner: the strong winds associated with tropical cyclones are represented by a turbulent shear flow over a flat surface, with a fluid of same density and viscosity as air. The sea spray droplets are modeled as small solid particles, which are dispersed in the flow, and can exchange momentum (but not heat) with the surrounding air. In this chapter, the governing equations for our simplified problem will be presented, along with the main assumptions defining our study framework.

Dispersed phase
Originally developed by Basset (1888), the solution for a sphere moving in a viscous fluid has been revisited numerous times since then, and improved by several authors who explored increasingly complex situations. Maxey and Riley (1982) suggested an equation of motion for a small sphere in a non-uniform flow.
This equation is at the base of the Lagrangian approach to deal with particles dispersed in a turbulent flow, and is given here in its original form: The indice j=1,2,3 corresponds to the three coordinate directions x 1 , x 2 , and x 3 , the mass of the sphere is defined by m p , and m f is the mass of fluid displaced by the sphere. The gravitational acceleration is represented by g j . Instead of the fluid velocity, u j represents the undisturbed flow field in this equation and v j is the particle instantaneous velocity, the material derivative d/dt actually following the object along its trajectory. By contrast, the derivative D/Dt corresponds to the material derivative following a fluid element (Maxey and Riley,1982), (Michaeledes, 2006). The acceleration of the sphere is given by the sum of the following forces: • the gravity net effect, • the fluid acceleration×m f , • the added mass, • the Stokes drag, • and the Basset history force.
The Faxen terms (i.e. ∝ d 2 p ∇ 2 u j ) are also included in the equation to take into account the curvature of the velocity profile around the rigid sphere. More recently, Ferrante and Elghobashi (2003) presented an up-dated version of the particle equation, adapted to the problem of multiphase turbulent flows where a large number of small spheres are suspended in a fluid. The equation they present can be applied to each particle to determine their respective acceleration: The first five terms on the right-hand side of the equation are easily identified as the forces presented in the original equation of Maxey and Riley. However the authors neglected the Faxen correction terms and introduced here the Saffman's lift force, which is due to the shear existing in the flow. The lift coefficient is defined as C L , L j is the direction cosine, and W is the magnitude of the relative velocity for a single particle.

Forces description
Gravity is generally the dominant force in presence, together with the Stokes drag. While the gravitational effects are not always considered in a large number of numerical studies, it has been observed that when taken into account in the equation of motion, the gravity would induce an anisotropy in the momentum transfer between the particles and the turbulence (Ferrante and Elghobashi, 2003), in the direction of the gravity vector g. On the other hand, the gravity tends to reduce lateral dispersion in multiphase flows (Sirignano, 2010) through a mechanism where the particles trajectories would not follow the eddies but cross them along their path (Crowe, 2006), (Ferrante and Elghobashi, 2003). Gravity also accentuates the preferential accumulation process regarding the dispersed phase being found mainly in the regions of low vorticity (Wang and Maxey, 1993), when compared to the case where its effect is neglected in the equation of motion.
Regarding the fluid acceleration, this term contains the fluid stresses from both the pressure and the viscosity of the fluid (Maxey and Riley, 1982). It can be understood as a net force applied to the sphere by the fluid surrounding it.
The added mass force, also called "virtual-mass" force, is due to the acceleration of the fluid surrounding a small sphere immersed in a non-uniform turbulent flow.
The volume of fluid being accelerated due to the sphere's motion is then equal to half the volume of the sphere, hence the term 1/2m f in front of the relative acceleration. It should be noted that the expression of the added mass force given here has been simplified by assuming that the particle Reyolds number is very small. When the local Reynolds number is finite, the original form of the added mass becomes: where the difference between the material derivatives is considered in the computation of the force.
The Stokes drag, fourth term found on the right-hand side of the motion equation, is the most important force in absence of gravity. It always acts in the direction of the flow (longitudinal force), and is proportional to the relative velocity between the sphere (or more generally the dispersed particles) and the carrier fluid. Its expression can be derived from the traditional drag force applied on a sphere: where C D is the drag coefficient, and u j is now the instantaneous fluid velocity at the particle location x p . For very small particle Reynolds number, the Stokes drag law can be applied and the expression of the drag coefficient becomes: For a sphere, the cross section A reads: leading to the drag force exerted in the x j direction by the sphere or the individual particles on the fluid: at their respective position x pj . This hydrodynamic force holds usually the main part of the interaction occuring in the two-way coupling between the particles and the turbulent flow.
The last term of the force balance, the Basset history corresponds to the diffusion of the vorticity around the moving sphere (in the case described by Maxey and Riley, 1982), which decays at a rate proportional to t −1/2 (Michaeledes, 2006).
Another interpretation would be that this term acts as a correction to take into account the transient character of the velocity field, in the case of creeping flows.
It translates as an additional resistance to the flow, though its effects decrease for either heavy and/or large particles (i.e. large Stokes number) (Crowe, 2006).
The Saffman effect can be explained as the force applied to a small sphere by the surrounding fluid due to the gradient of the streamwise fluid velocity in the wall-normal direction (shear flow). The lift force is generally small compared to the other terms and can be neglected in most cases. More specifically, it is much weaker than the conventional drag force acting in the streamwise direction of the flow.
Nevertheless, in the specific case of turbulent channel flows where the turbulence is bounded by parallel planes, the Saffman lift force has been shown to impact the deposition rate of the particles near the walls. In their paper, Marchioli and Soldati (2002) observed that the particle deposition at the walls is enhanced when the particles move faster than the fluid, and reduced when they are slower than the carrier fluid. Moreover, the lift's effects are weaker when the particle inertia is important (i.e. large and/or dense particles), and remain confined to the viscous sublayer of the flow. In the case where the lift force is omitted to model the dispersed phase equations, the particle fluxes towards the wall have been seen to decrease slightly (Marchioli and Soldati, 2002), without modifying qualitatively the results obtained. It should be noted that the treatment of the Saffman's effect in a turbulent channel flow is not straight-forward, and requires particular attention when dealing with the particle-wall interactions, especially at higher particle Reynolds numbers.

Assumptions
In regards to the work of Maxey and Riley (1982), we define the characteristic length scales, respectively for the undisturbed flow and for the small sphere as L and d p (the diameter of the sphere). The corresponding characteristic velocity scales are U for the flow, and W for the relative velocity between the sphere and the fluid around it. We make three major assumptions regarding their equation of motion: • the moving sphere should be very small compared to the length scale of the flow: d p /L 1 (pointwise approach), • the local (particle) Reynolds number should be very low: d p W/ν 1, • and the velocity gradients of the mean flow should be small as well: (d 2 p /ν)(U/L) 1 (shear Reynolds number condition).
In our study however, where droplets of seawater are suspended in the air, supplemental hypotheses can be made to simplify the equation of motion for the dispersed phase further. More specifically, the following assumptions are made: • the particles are rigid, i.e. no deformation is allowed, and the shape and dimensions of the droplets remain constant, • the particles are spherical, i.e. for simplicity, no shape effects are considered, • the particles are heavy, i.e. the density ratio between the water "particles" and the air is large (ρ p /ρ f = 1000).
In the framework of DNS of turbulence, the shortest length scale of the flow is usually the Kolmogorow length scale η K . The original constraint on the size of the particles becomes then d p /η K 1, and the condition of low particle Reynolds number, redefined as Re p = d p | u j − v j | /ν, still holds so that the Stokes law applies.

Scaling
In most practical cases, not all the forces described in the previous section need to be accounted for to determine the particle velocity and to track the position According to the third assumption suggested by Maxey and Riley (1982), (d 2 p U/νL) is small. In our case, the two velocities U and W are at most of the same order of magnitude, so that their ratio tends to 1. Therefore, the fluid acceleration is always negligible. The Basset history term, when compared to the Stokes drag, is of order O(d 2 p U/(νL)) 1/2 , while the added mass term compared to the drag force leads to terms of order O(d 2 p U/(νL)). Again looking at the third assumption, we can conclude that the Stokes drag dominates both the added mass force and the Basset history in the motion equation. The lift force when scaled by the drag force is of order of the particle Reynolds number, which we assumed to be small: we can thus ignore the lift effects on the particles in the present study. The Faxen terms, appearing in the equation of Maxey and Riley (1982), are second-order terms and can be neglected safely under the assumption that the particles are much smaller than the characteristic flow length scale. Michaeledes (2006) presents another way to scale the history term and the added mass relative to the drag by expressing the particle equation in a dimensionless form. There, the velocities are scaled by the characteristic velocity of the flow and the time scale for the sphere τ p , introduced in the next section. The author shows that in that case, the Basset history scales as (ρ f /ρ p ) 1/2 , while the added mass scales as ρ f /ρ p . By noting that the particles are much heavier than the fluid they are suspended in, these two forces can again be safely neglected in front of the Stokes drag in our context.

Reduced equation
Following the scale analysis, we consider only the drag contribution in the force balance, and we ignore the contribution of the fluid acceleration, the virtualmass, the Basset history, and the Saffman lift effects respectively, as well as the Faxen terms. The simplified version of the particle equation of motion can then be applied to our case: In the absence of gravity, the drag balances exactly the particle acceleration, and we can write: with m p being the product of the particle density and the volume of one spherical particle, m p = ρ p πd 3 p /6. We note F pj as the force acting on the particle p in the x j direction and by the fluid. The particle equation of motion can then be expressed again, this time by including a specific time scale, named particle response time 15 The particle response time is a crucial characteristic time scale of multiphase flows which represents physically the time necessary for the dispersed phase to respond to the fluctuations of the flow velocity. The larger τ p is, the more inertia the particles will display. At the limit where τ p is very large, the inertia of the particles is so large that the particles behave like ballistic projectiles, while when τ p is close to zero, the particles act as passive tracers in the flow, because they are able to follow the fluid freely everywhere in the domain. Furthermore, this time scale is crucial to determine the behaviour of the dispersed phase, when comparing τ p to a typical time scale of the turbulent flow.

Carrier phase
Our current model to describe the fluid governing equations is based on the work of Ferrante and Elghobashi (2003), which dealt with the physical mechanisms of two-way coupling in particle-laden isotropic turbulence. In a particle-laden incompressible flow, the Navier-Stokes equations can be applied for the fluid phase: along with the continuity equation: where u j are the fluid velocity components, and ∂p/∂x j is the pressure gradient.
The density of the fluid is denoted ρ f , and µ is the dynamic viscosity of the fluid.
The two-way coupling between the phases is related to the term −F j , which is the net force exerted in the x j direction by M particles (that exist in a unit volume of the fluid), and is computed from F pj being the steady-state drag force described in the previous section for a single particle.

Numerical Methods
In this chapter, we will present in detail the methods used to simulate the interactions between small inertial particles (standing for the sea spray droplets) and the turbulence of a channel flow. Our framework lies on an Eulerian-Lagrangian approach, where the fluid's motion is solved in the Eulerian frame through the Lattice Boltzmann Method (LBM) while the dispersed phase motion equations introduced in the previous chapter is solved in a Lagrangian way via a concise numerical scheme. This scheme is then implemented together with the core of the LB code, in order to solve simultaneously the particle equations and the Navier-Stokes equations at each time step of the DNS runs. In addition to the description of the numerical approach, the geometry of the problem as well as the important non-dimensional parameters of our study will also be discussed here.

Lagrangian method
As explained in the previous chapter the particle equation, or motion equation, relates the particle instantaneous acceleration to the sum of the forces applied on the same particle. Through scaling arguments we have shown that in absence of gravity, only the drag force is important to compute the particle velocity. The equation can be solved at each time step, and integrated in time to obtain the particle position x p (t): To determine the particle velocity, however, the knowledge of the fluid velocity at the particle position u j [x pj (t), t] is essential. The fluid velocities are computed via the LBM in the Eulerian framework, and are only known at the grid nodes of our computational domain. Therefore, to derive the values of the fluid velocity field at each particle position, an interpolation of the velocity field at each time step is required to retrieve the drag force (i.e. the coupling term between the carrier fluid and the dispersed phase F j ), and track its evolution over time.
Over the course of this project, we have tried several interpolation schemes with various degrees of accuracy to calculate the fluid velocity field at the particle position. We have finally chosen the method proposed by Lekien and Marsden (2005), which uses a tricubic interpolation in three dimensions. This local interpolation method is based on the determination of a 64x64 matrix that relates the derivatives at the corners of an element to the coefficients of the tricubic interpolant for this element, and presents two main advantages: on one hand it uses only the neighborhing points of the element instead of the whole dataset to determine the fluid velocity; on the other hand, a unique set of coefficients for the velocity interpolation is determined once and stored for further usage, which saves both time and computational ressources. We will briefly explain their method, but more details are described in the original paper.
For simplicity, we use the same notation as the original paper of Lekien and Marsden (2005). Let us consider a small cubic element as part of the computational mesh. The figure 2 illustrates the element considered in this approach, from the original paper. In our case, we take a cubic element of width equal to one grid node to avoid unnecessary complications. The fluid velocities, as well as their first derivatives in space, their second order mixed derivatives and their third order mixed derivative are computed at each corner of the element, such that we have a set of eight values, computed at the eight points denoted p i , i = 1, 8; and f being the velocity function here, i.e. u j (x, y, z, t): 19 The function f can be expressed as a polynomial of the form: For the sake of computation, the 64 coefficients a ijk of f can be expressed under the vectorial form through the transformation: The concept here is such that p 1 corresponds to the origin of the cube (0,0,0), and the coordinates of the points p 2 , and p 8 are respectively (1,0,0) and (1,1,1).
Figure 2: Schema of an element for interpolation in three dimensions (Lekien and Marsden, 2005) By computing f and its derivatives at the eight points (0, 0, 0), (1, 0, 0), ..., (1, 1, 1) , we thus obtain a vector containing 64 variables, which we will call b i . Then, using the expression of f introduced above, it is possible to relate the vector b i to the vector α i through a system of linear equations, such as: Here, B −1 is the 64x64 matrix mentioned above, which was preliminary determined, and implemented in the GPU code to interpolate the fluid velocity field within the computational domain. The main advantage to use this method lies on the fact that the determination of the polynomial coefficients happens once, and does not need to be actualized at each time step. To our knowledge, this method has not been used before for similar particle-laden turbulent flows problems. Most previous studies using more traditional numerical methods, such as the pseudo-spectral method to solve the carrier phase equations (Soldati et al., 2002(Soldati et al., -2012, , have relied on a classic 6th order Lagrangian polynomial interpolation scheme to derive the fluid velocities at the particle positions.
Once the fluid velocities are known, the drag force can be easily deduced from the expression developed in the previous chapter: This force is computed for each particle suspended in the fluid, and is up-dated at each time step with the fluctuacting velocities. We recall that it corresponds to the interaction term, or feedback from the dispersed phase to the carrier phase, and that it is the key to the two-way coupling approach we chose where the phases are expected to exchange momentum with each other. The equation 2 can be written again, bearing in mind that it corresponds to an application of Newton's second law: with a j = dv j /dt and the right-hand side term being equal to [u j (x pj ) − v j ]/τ p .

21
The particle instantaneous acceleration can therefore be known readily: From this, the particle velocity can be computed at the next time step through a first-order Euler-forward scheme: The up-dated particle position is obtained via a time integration taking into account both the velocity and the acceleration of each particle: It should be noted that the time steps used in the LBM are very small, which allowed us to compute both the particle velocity and the particle acceleration with reasonable accuracy. The last step of the particle-tracking scheme consists in smoothing the drag force in the domain, so that the fluid governing equations (14) can be solved at the next time step while taking into account the feedback contribution from the particles. For stability reasons, the drag needs to be extrapolated after being computed in a pointwise manner at the particle positions. Over the course of this study, we tried several extrapolation methods: we first used weights proportional to the distance of a particle to its eight neighboring grid nodes, then looked at a radial smoothing scheme for the feedback force over a distance proportional to the magnitude of the force. However we chose to use a simpler trilinear extrapolation, with constant weights equal to 1/8. This basic scheme was preferred to the more elaborate ones due to the large number of inertial particles dispersed in the flow: it allowed us to save computational time by imposing constant weights instead of deriving their values for each particle, at each time step. This spatial redistribution of the force prevents peak values of the drag from occurring across the computational domain.

22
To summarize, we distinguish five stages in the Lagrangian tracking method: • the interpolation of the fluid velocities, • the computation of the drag force, • the determination of the particle velocities, • the time integration to obtain the particle positions, • the extrapolation of the drag.
These steps are implemented in an autonomous code to simulate the advection of the particles suspended in the turbulent flow, and are repeated at each time step.
This code communicates to the other parts of the GPU code by using the fluid velocities as an input and giving the smoothed drag over the domain as an output to be applied in the governing equations solved by the Lattice Boltzmann method (LBM).

Lattice Boltzmann Method
The Lattice Boltzmann method takes its roots in the kinetic theory of gases, and has evolved from the Lattice Gas Automata models developed in the 1980s to become an enticing alternative to traditional numerical methods to solve various flows and engineering problems. In this section, we will describe briefly in what consists the LBM. More details can be found in the reference book by Succi (2001) on the background of the Lattice Boltzmann equation, and its applications in computational fluid dynamics.
According to Succi (2001), we can define the kinetic theory as "the branch of statistical physics dealing with the dynamics of non-equilibrium processes and their relaxation to thermodynamic equilibrium". At the microscopic level, let us first imagine an ensemble of fictitious fluid particles. These particles can be considered as pointwise, and respect the principle of continuum: they contain a large number of molecules, but are still small compared to the macroscopic scale, such that the fluid in its globality can be seen as a continuum medium. We introduce now the distribution function f (x, p, t) as the probability of finding a molecule around the position x at a time t with the momemtum p. Boltzmann (1872) developed an equation to describe the evolution of this distribution function, which is known now as the Boltzmann equation (BE): where Ω represents the collision operator, F represents the sum of the external forces applied to the particles, and m is the mass of the fluid particles. The lefthand side of the equation corresponds to the advection of the function f , and is also called the streaming part of the BE.
The earlier numerical models solving the dynamics of these particles were the Lattice Gas Automata, where time, space and the particle velocities were all discrete. The fictitious fluid particles were found on the nodes of the lattice and at each time step, they could either move to the nearest node in the direction of their velocity, or collide with the neighboring particle, since two particles could not occupy the same node at the same time. Hence, the two parts of the BE correspond to the two possible configurations for the fluid particles (propagation and collision), as shown in the figure 3.
In the original LGA models, the distribution function was defined as a boolean variable. However, in the Lattice Boltzmann method the distribution function is actually the average of this boolean variable. Furthermore, to reduce computational costs, particles are only allowed to move along certain directions and with given speeds, so that the set of discrete particle velocities becomes finite. The common BE (in absence of external forces this time) can then be expressed as: with i = 0, 1, ...Q, Q being the number of directions for the particle velocities, e i the particle velocity and Ω i still the collision operator. A finite difference discretization in time and space leads to the Lattice BE: with ∆t being the time step. This equation applies under the condition that the particle speed of propagation is set to c = ∆x/∆t = 1, with ∆x being the grid spacing of the lattice (Janssen, 2010). To determine exactly which set of particle velocities are being used in our model, we refer to the classic convention used in the Lattice Boltzmann community where models are named according to their dimension and the number of discrete velocities involved. For instance, a vastly used model is the D2Q9 model, as a two-dimensional LB model with a set of nine discrete velocities. In our case, since our problem is 3D, we use the D3Q19 model where the total number of particle velocities allowed in the velocity space are 19. The figure 4 illustrates the repartition of the velocities in space for the D2Q9 and the D3Q19 models. One of the remarkable feature of the LBE Figure 3: Schema of the propagation and the collision steps (Janssen, 2010) is that the incompressible Navier-Stokes equations can be recovered through the Chapman-Enskog expansion (1916)(1917). Indeed, by expanding the distribution function f i , it can be demonstrated that the microscopic LBE converges towards the macroscopic Navier-Stokes equations, such as: and where M a is the Mach number (discussed in the next section) and Kn corresponds to the Knudsen number, defined as the ratio between the molecular mean free path and the shortest length scale of the flow (Succi, 2001).
Figure 4: Schema of the D2Q9 and D3Q19 (Janssen, 2010) For simplicity the external forces (gravity and particle feedback) have been omitted in this version of the Navier-Stokes equations. It should be noted that in contrary to traditional pseudo-spectral methods, the pressure is solved via an equation of state, rather than through the computational-demanding Poisson equation.
Macroscopic variables, such as the fluid density and the momentum can be also retrieved from their microscopic counterparts as follows: The conservation of mass and momentum in the computational domain can be translated as the following conditions: and Ω i e i = 0 respectively.

The collision operator contains information about the interactions between fluid
particles. In particular in the BE, it focuses on the localized binary collisions occurring between particles by assuming that the fluid is dilute enough to avoid more complex interactions: only short-range binary collisions are considered in this framework. In the LBM, there are several possible expressions for the collision operator Ω, and the choice of one formulation over another depends usually on the problem considered. Without going into further details, we will just specify that we use a multi-relaxation time (MRT) model to represent the collision operator in our problem. It is known to have a better stability and to provide more accurate boundary conditions than the standard BGK (for Bhatnagar, Gross and Krook) model by involving higher moments of the distribution function in the expression of Ω i .

Geometry
Within the framework of our study, the choice of the flow configuration was made so that we would be able to compare our results against several published 27 studies. A 3D Poiseuille turbulent flow, bounded by two flat planes, was simulated with and without small particles, to observe the interactions between the dispersed phase and the fluid. For simplicity, the channel planes are smooth and no-slip boundary conditions are imposed at each plane. In the LBM, no-slip boundary conditions are implemented under the form of a "bounce-back" scheme: the fluid particles bounce off the wall after reaching the boundary, in the same direction that they came from (Janssen, 2010), and while conserving their mass and momentum exactly. More specifically in our simulations, a second-order bounce-back scheme with a higher accuracy in space was used for the carrier phase: by placing the channel planes half-way between two lattice nodes, the last row of nodes in the vertical direction taken into account for our computations were at a distance d = 0.5∆z, and a row of "ghost nodes" located beyond the walls was then considered (Janssen, 2010). The flow is periodic in the homogeneous directions, i.e. in the streamwise and the spanwise directions of the flow. Regarding the inertial solid particles, periodic boundary conditions were also implemented: once a particle leaves the domain in one of the horizontal directions, it is reinjected on the other side of the domain, with the same velocity.
The flow is driven by a mean pressure gradient in the streamwise -x direction, which is implemented as a body force in the LBM code. The figure below schematizes the geometry of our problem. The dimensions of the planes were chosen so that the flow can be considered homogeneous in the streamwise and the spanwise directions (Bespalko, 2011). The main constraint in the choice of the domain size lies in the fact that the largest structures of the turbulence should be contained within its whole geometry. Previous published results by Kim et al. (1987) showed that a size of the computational domain of 8δ x 4δ x 2δ was sufficiently large to ensure that the two-point correlations for the velocity and the pressure fluctua-28 tions decay to zero. They also showed that grid spacing should be smaller than δ/90 to accurately resolve the viscous sublayer. Then, the optimal dimensions of the numerical domain would be 720 x 360 x 180. It should be noted that this configuration requires at least two GPGPUs to ensure enough memory is available to save all the data (i.e. the fluid and particle velocities, as well as the particle parameters). As of now, the computational domain limited to 260x260x180 nodes, in the streamwise, spanwise and wall-normal directions, respectively. Therefore, our results may be slightly affected by the fact that the largest turbulent eddies are not well resolved.

Dimensionless numbers
In this section, we will introduce the parameters that determine the framework of our simulations. These non-dimensional numbers are important to make sure that our simulations correspond to some real-life cases by matching their values respectively between the physical space and the computational domain. Besides the Mach number which is intrinsic to the Lattice Boltzmann method, there are four main parameters to consider when dealing with particle-laden turbulent flows.
A fifth one can be counted in as well if the gravitational force is taken into considereation in the problem.  Volume fraction The volume fraction, denoted φ v , allows to determine how much of the dispersed phase is present in the flow. It is the ratio of the volume occupied by the particles to the total volume (fluid+particles):

Mach Number
Alternatively, the mass fraction can be used as well to specify the global concentration of particles in the flow. In a similar way as the volume fraction, it then corresponds to the ratio of the mass of all the particles dispersed to the total mass of the system: Since the focus of our work is the possible interactions occuring between the air turbulence and the sea spray droplets at the ocean surface, our interest is confined to the two-way coupling regime where the carrier phase (air) and the dispersed phase (spray) can exchange momentum with each other (Elghobashi, 1994). This implies that the range of concentrations we will investigate is 10 −6 for the lower bound (below this value, the suspension is too diluted to allow any effect of the droplets on the turbulence) and 10 −3 for the upper bound (above this value, the suspension is too concentrated, and the droplets start to interact with each other: coalescence, breaking up...).
Density ratio As its name indicates, the density ratio is expressed as the ratio of the density of the dispersed phase over the density of the carrier phase: When considering the forces that the carrier fluid apply to the dispersed phase, the density ratio distinguishes between the forces that can be neglected, and the ones that are relevant to the problem (Maxey and Riley, 1982). For instance the particle inertia, characterized by the particle response time, becomes predominant when D is very large. In our case, the density ratio is equal to 1000, which means that the particles can be qualified as heavy. Consequently, the Stokes drag is much larger than the other forces like the added mass and the Basset history term. D will remain unchanged throughout our study, unless specified.

Results
The impact of the sea spray droplets on the near surface turbulence is

Simulation parameters
The definition of the simulation parameters is thoroughly discussed in the PhD thesis of Bespalko (2011). Here, we present the complete set of parameters necessary to perform our DNS runs.
In contrast with free shear turbulence, wall-bounded turbulence exhibits two main regions: a region close to the wall where viscous effects are important, and another region further away where the turbulent stresses dominate the viscous stresses. The two regions are called the viscous wall region and the outer layer respectively (Pope, 2000). The important notion to retain is that the turbulence is governed by different scales depending on the regions one is interested in. For instance in the viscous layer, the key parameters are the kinematic viscosity of the fluid ν, the density of the fluid ρ f (or conversely the dynamic viscosity of the fluid µ which depends on both the density and the kinematic viscosity), and the wall shear stress τ w . From these three quantities, it is possible to derive characteristic length, time and velocity scales, as well as a non-dimensional number characterizing the turbulence state at the wall (i.e. the friction Reynolds number from the previous chapter). The velocity scale was briefly introduced before: it is the friction velocity, here formally experessed as: The viscous time scale was also presented when dealing with the viscous Stokes number, and is defined as τ v = ν/u 2 τ . Lastly the viscous length scale is: According to KMM (1987), one of the requirements of a turbulence DNS is that the grid resolution should be fine enough to resolve the smallest length scales of the turbulence, and in particular it should be of the order of the Kolmogorov scale η K . The Kolmogorov length scale is estimated to be equal to about 2 wall 35 units, when normalized by the viscous length scale δ v : η + K = η K /δ v ≈ 2. (Here, the superscript + denotes a quantity normalized by the viscous layer scale.) In channel flows, the grid resolution should be at least equal to the normalized Kolmogorov length scale in the vertical (wall normal) direction. We therefore impose the normalized grid resolution to be ∆z + = η + K = 2 in our configuration.
The grid resolution in the other directions need not be as small as the Kolmogorov length scale. However, the LBM is normally designed so that the lattice where the time, space and velocities are discretized is uniform in all directions. This means that the lattice is cubic, with ∆x + = ∆y + = ∆z + . As a result, the flow is over resolved in the horizontal directions, while being appropriately resolved in the wall-normal direction.
In the Lattice Boltzmann Method the length, time, and mass scales are defined differently from the ones in the physical space. Here, the symbols L, T and M will be used to characterize the dimensions of length, time and mass respectively in the lattice system. In particular, the grid spacing is set 1 in the lattice length scale. Hence, the viscous length scale in the lattice unit system is equal to 0.5L.
To determine the value of the friction velocity in the lattice space, it is necessary to keep in mind that the LBM is a compressible numerical method on which we imposed a given value for the Mach number so that the problem we are interested in remains in the incompressible regime. In the previous chapter, we set up the Mach number to be M a = 0.2. By definition, the speed of sound in the D3Q19 model is given as (Bespalko, 2011;Janssen, 2010): which leads to a centerline mean velocity U 0 ≈ 0.115L.T −1 . According to the friction law, derived by Bespalko (2011): The von Kármán constant κ, and A are respectively equal to 0.4 and 5.5. In the expression above, the mean velocity at the centerline of the channel has been normalized by the friction velocity u τ , and Re τ is the friction Reynolds number defined in previous sections. With all the parameters known, the value of the friction velocity can be deduced from the friction law: u τ ≈ 6.25x10 −3 L.T −1 . Going back to the expression of the viscous length scale as a function of the kinematic viscosity and the friction velocity, the fluid viscosity is then: ν ≈ 3.12x10 −3 L 2 .T −1 .
Once the friction velocity is known, the wall shear stress can be calculated from the expression given earlier in this chapter: The determination of the wall shear stress is necessary to compute the mean pressure gradient required to sustain the turbulent Poiseuille flow along the streamwise direction. In the configuration of channel flows, the mean pressure gradient should be exactly balanced by the sum of the stresses, at every point of the computational domain. In particular, at the wall where the Reynolds stresses are zero, the only component of the stress is the wall shear stress, and we can thus write: Here, the value of the mid-channel height in the lattice space is obtained from the definition of the friction Reynolds number: As specified by Bespalko (2011), the pressure gradient is implemented in the LB code as a body force (−∂p/∂x = ρ g), instead of a true pressure gradient to avoid dealing with mean pressure gradients in the flow.
In the following section, we will present the results for the clean channel case, compared to the benchmark published by KMM in 1987. The main statistics will be shown here, as well as the results for different domain sizes.

Turbulent Poiseuille flow: Validation 4.2.1 Two-point correlations
In order to model the interactions between small particles and turbulence, we need to ensure first that the implementation of turbulence is correct before introducing particles in the flow. Our original intent was thus to reproduce the results published by the reference article of KMM (1987) on turbulence statistics in a fully developed channel flow at low Reynolds number. In their paper, the authors observed that a computational domain of dimensions 8δx4δx2δ was sufficiently large to ensure that the largest structures of the turbulence were contained within the whole geometry by examining the two-point correlations for the velocity and the pressure fluctuations. The two-point correlations constitute a simple way to obtain information on the spatial structure of a random field (Pope, 38 2000). Kim et al. (1987) showed that the correlations were indeed decaying to zero for large separations, supporting their choice of the domain size.
Regarding the LB configuration where the grid has to be isotropic, this would imply that the optimal dimensions for the numerical domain would be 720x360x180 nodes, in the streamwise, spanwise and in the wall-normal directions, respectively.
However, it should be noted that this configuration would require at least three GPGPUs to ensure that enough memory is available to save all the data (i.e. the fluid and particle velocities, as well as the particle parameters) for each run. In this project, we were able to use only one GPU to run our DNS, which limited drastically our computational resources. We maintained the constraint on the grid resolution in the wall-normal direction, so that our results could be compared to the study of KMM. However, our domain was cut shorter in the x-and -y directions: the geometry used in our work was 260x260x180 nodes instead of 720x360x180 nodes. Knowing that the channel we considered for this study was smaller than the optimal size prescribed by KMM, we computed the two-point correlations for the velocity fluctuations to observe the discrepancies due to the smaller geometry.
The general expression to compute the two-point correlation is given by Pope (2000): where r is the separation between the points, and x is the position vector. It is common to view turbulent flows as a superimposition of a mean flow and some fluctuations over time and space. The Reynolds decomposition allows one to distinguish between the mean and the fluctuating components of the flow variables, such as the velocity and the pressure. We can then write:

39
The brackets designate an ensemble average, performed as follows: where the sum of the velocities is calculated over time and across the horizontal plane at the height z of the computational domain. n t is the number of snapshots considered for the time average, while n x and n y correspond to the number of grid points in the x-and y-directions respectively. The fluctuating part, denoted with a prime, can be deduced after computing the average in a straight-forward manner: The figure 6 shows the two-point correlations for the fluctuating velocities at different heights in the channel. Each time, R ij was normalized by the average velocity at those respective heights. Our results can be compared with the results in the reference paper, shown in figure 7.
Since our domain is only 36% of the optimal length in the streamwise direction, a comparison would be relevant only up to around x/δ ≈ 1.4 in the reference paper. Even with the largest separation in our computational domain, the correlations have not reached zero yet. Nevertheless, the overall pattern is quite similar between the two results. Our R uu is slightly higher than the original results at both channel heights. But, our R vv and R ww show good agreement with the reference results: R vv reaches zero near the wall, and is slightly negative near the channel center, while R ww is positive near the wall and negative but reaching towards zero at the mid-channel height. Despite our computational limitations with a smaller domain on a single GPU, it appears that the two-point correlations agree reasonably well with the results of KMM's paper, up to the point of the largest separation δx in our problem. Therefore, our results are expected to converge towards the benchmark study with a larger computational domain.

Mean velocity and Reynolds stresses
In order to investigate the convergence towards the results published in the reference paper, we ran two DNS cases, with varying gridsizes. In case A, the domain size corresponded to the one used in this whole project (which is about the largest size reachable on one GPU alone with particles), namely 260x260x180 nodes, while the channel was larger in case B, with dimensions of 350x300x180 nodes (which is about the largest size reachable on one GPU alone without particles). The figure 8 shows the mean streamwise velocity, as well as the normalized Reynolds stresses, defined as u i u j /u 2 τ (Pope, 2000) for the cases A, B, and KMM. The mean streamwise velocity, normalized by the friction velocity is plotted in a semi-logarithmic scale as a function of the distance to the wall, which has been normalized with the viscous length scale. The thin dash lines correspond to the viscous law u + = z + , and the Nikuradze log law (Pan and Banerjee, 1996) u + = 2.5 lnz + +5.5 respectively. We observe a very good agreement between the reference paper and our results for the streamwise mean velocity. Concerning the stresses, the largest discrepancies are seen for the streamwise terms. As expected, the results in the wall-normal direction are the closest to the KMM's results among the three directions since we respected the resolution constraint in the wall-normal direction.
Generally, the larger domain results show a better agreement to the benchmark study than the smaller domain results, which corroborates the assumption that our simulations would converge towards the published results of KMM with a larger domain size.

Higher-order statistics
In addition to the mean streamwise velocity and the Reynolds stresses, we have examined the higher-order moments of the velocity fluctuations: the skewness (third order) and the flatness -also called kurtosis (fourth order). The figure 9 illustrates the skewness and the flatness plotted along the channel height. The expressions for both moments are given here: Again, our results are in good agreement with the benchmark results of KMM (not shown here). From the two-point correlations to the fourth order moment, our results are in good agreement with the reference paper by Kim et al. (1987), even though our computational domain is a truncated version of the optimal channel described in KMM. This also suggests that the LBM can reproduce accurately DNS of turbulent Poiseuille flow, provided that the grid resolution and the gridsize requirements are met. Since the grid is by default isotropic in LBM, the flow is over resolved in the wall parallel directions, which indicates that the LBM is not the most numerically economical method to solve turbulent wall-bounded flows. Nevertheless, in regards to turbulent particle-laden flows, the efficiency of the LBM can be appreciated in terms of computational time and memory required to deal with such flows.
The particle equation of motion is solved simultaneously with the LBM at every time step, and the small particles are tracked in their trajectory over time, while exchanges of momentum are allowed between the dispersed phase and the carrier fluid. In the next section, the low-order statistics of turbulence will be presented in the case where particles are introduced in the channel flow. Additional results will also be discussed regarding the modulation of turbulence by inertial particles.

Turbulent particle-laden flows
We recall that there are four non-dimensional numbers that govern the dynamics of turbulent particle-laden flows: the friction Reynolds number, the Stokes number, the volume/mass fraction and the density ratio. In this study, we kept the Reynolds number and the density ratio constant, and varied the concentration and the size of the particles (i.e. their Stokes number).  The main numerical results of this study will be presented in the following sections of this chapter, with the following convention: most figures will be plotted as series of curves, obtained first by varying the mass fraction (at constant Stokes number), and then by varying the Stokes number (keeping the mass fraction constant). The legends, which will be consistent for the whole chapter are described in table 2.  The disparity with the clean channel becomes more pronounced when a larger number of particles is dispersed in the flow.
Regarding the cases I, III, V and VI where the Stokes number is progressively increased from 10 to 720, the streamwise velocity fluctuations are higher than for the clean channel, though the difference is negligible for the case VI.
The enhancement of the streamwise fluctuations is particularly evident for the cases III and V. Dampening of the r.m.s in the y-and z-directions can be seen in all cases. Moreover, the case I with the smallest particles displays an interesting behavior: the curves for v rms and w rms start at the wall between the cases III and V, but end up crossing the curves after the viscous wall region.
Overall, we may conclude that there is a global effect of φ m and of St + on the lower-order statistics of turbulence, though the underlying mechanisms of turbulence modulation may be different whether St + or φ m are changed in the DNS runs.
According to Zaho et al. (2013), it is thought that the inertia of the particles play an important role in modifying the low-order moments of turbulence, at any

Force balance
It is important to make sure that the two-way coupling is correctly implemented: the feedback force, which is included in the Navier-Stokes equations from the previous chapters, exchanges momentum between the fluid phase and the dispersed phase. We verify here that the conservation of momentum is satisfied, and focus on the particles' role in the global force balance. The figure 13 shows the repartition of the various forces in presence in the channel flow as a function of the normalized distance from the wall. In the case of the clean channel, the force balance consists of the viscous force and the Reynolds force compensating the mean pressure gradient exactly: Note that while many studies examine the stress balance and compute the viscous and Reynolds stresses in turbulent Poiseuille DNS, examining the force balance is more relevant in the case where particles are dispersed in the computational domain. We emphasize that both formulations are equivalent, and we can relate the forces described in the equation above to the viscous stress and the Reynolds stress as: In the case of particle-laden flows, the force balance becomes: The results of the force balance are shown in figure 13. Note that the legend in this figure differs from the previous plots: the solid lines in all cases correspond 51 to the viscous force -regardless of the color, the dash lines represent the Reynolds force, and the dotted lines show the particle terms, taken in the streamwise direction. The red starred line is the sum of all forces in presence, and the black thin horizontal line found at -4.34x10 −7 M.L −2 .T −2 is the mean pressure gradient.
The first plot corresponds to the clean channel case, where we can see that the sum of the forces exactly compensates the mean pressure gradient applied to the flow, proving that there was no sink or source of momentum in the channel flow. The second plot illustrates the evolution of the repartition of the forces among the viscous force, the Reynolds force and the particle feedback force for the cases II-blue, III-green and IV-red (increasing concentration), while the third plot shows a similar family of curves, for the cases I-blue, III-green, V-red (increasing Stokes number). The last plot on the lower right corner shows the averaged force components, for the case with the largest concentration of particles (case IV).
As expected for the clean channel case, the absolute value of the viscous force is maximal in the viscous layer (Pope, 2000), and reaches zero beyond the buffer layer (z + ≈ 30). The Reynolds force has an opposite sign to the viscous force, and compensates it exactly by reaching the value of the mean pressure gradient at the point where the viscous effects vanish. The effect of the increasing mass fraction on the force balance can be seen in the second figure in the upper right corner: inside the viscous sublayer the magnitude of both the viscous force and the Reynolds force diminishes progressively, while the particle force F x increases at the same time for larger φ m . A small discrepancy between the sum of the forces and the mean pressure gradient is observed near the wall, which is due to the treatment of the no-slip boundary conditions in the LBM. Besides this numerical artifact, however, it is clear that the force balance is still satisfied in presence of particles.
On one hand, it should be noted that the Reynolds force decreases faster than the viscous force near the wall for increasing mass fraction. On the other hand, the maxima locations for the particle force for cases II, III and IV are found closer to the wall than the maxima of both F visc and F R , which indicates that the particles impact the turbulence mainly close to the wall, rather than in a homogeneous fashion over the whole domain. While its magnitude can reach up to five times the value of the mean pressure gradient, it remains very localized in the domain. Averaged particle force, case IV z + Figure 13: Force balance: in clean channel (upper left); varying mass fraction (upper right: case II-blue, III-green, IV-red); varying Stokes number (lower left: case I-blue, III-green, V-red); and averaged particle force components (lower right, case IV). Solid line: F x , dash line: F y , dash-dotted line: F z Outside the viscous sublayer the particle force is small but remains finite, and reduces the Reynolds force accordingly. The results with different particle sizes (Stokes number) is shown in the third panel. Again, the magnitude of the forces are smaller than for the clean channel case. However, the evolution of the forces is not monotonic: the curves for the case I are found between cases III and V (Reynolds force) and case V (viscous force). The smaller impact is seen for the case VI (not plotted here), which is expected since the less amount of particles is released in the channel flow.
From the figure 13, we may conclude that regardless of their size or concentration, small inertial particles allowed to interact directly with the carrier fluid impact the distribution of the forces by dampening the Reynolds force (across the channel) and the viscous force (within the viscous wall region).
Lastly, the fourth figure on the lower right corner illustrates the averaged particle force across the channel for Case IV. As expected, F y is always zero. It is interesting that the mean particle force in the wall normal direction F z is not zero near the wall, although its magnitude is much smaller than that of F x . Zaho et al. (2013) investigated the effects of the particles on the turbulent kinetic energy (TKE) budget, focusing on the energy transfer between the fluid and the solid phases, and on the dissipation in the wall turbulence. To take into account the particles into the TKE budget, we derive the TKE equation from the Navier-Stokes equations (eq.14) presented in the previous chapter. The TKE is defined as k = 1/2 u i u i , and the expression of the conservation of the turbulent kinetic energy is given here:

Production and dissipation: the TKE budget
Although the notation might differ from the expression given by Pope (2000), we can distinguish seven terms in the balance of TKE, namely: • the material derivative of the TKE, • the pressure transport, • the turbulent transport, • the viscous diffusion, • the production, • the dissipation rate, and • the particle production.
In the case of a statistically steady channel flow, the TKE is conserved and the material derivative of k, which corresponds to the sum of the local rate of change and the advection of k by the mean flow, is equal to zero. Since the gravity has been neglected in this problem, the buoyancy does not intervene in the TKE balance.
Here, we are interested in the change in the repartition of the production terms and the dissipation rate. We will note the TKE production P , the particle production P p and the dissipation rate ε. In wall-bounded turbulent flow, the production term P can be reduced to: The figure 14 shows the mean shear production and the particle production terms, as well as the dissipation rate, for the cases I to VI, with the clean channel results as baseline. The left-hand side of the figure represents the cases II, III and IV. We can observe here a reduction of the production and the presence of a non-zero term P p due to the feedback of the particles. The magnitude of the dissipation rate also decreases with increasing mass fraction. The top and bottom panels show that the particle production is small compared to the TKE production, and negative across the channel height for the cases III and IV. The fact that P p becomes negligible in the case II is due to the very low concentration of particles dispersed in the flow.
The right-hand side panel shows that the particles of various diameters impact in a more complex manner the production and the dissipation of TKE, although both the production and the dissipation are reduced by the particles. The impact is significantly reduced for the case VI (St + = 720, φ m = 0.25) since the number of particles becomes very small.
In the last figure, in the bottom right corner, the particle production changes sign: it starts negative for the cases I and III, then is positive across the channel in the case V, and finally is almost equal to zero for case VI where very few particles were dispersed in the flow (like in case II).
In summary, the introduction of small particles in the channel flow modifies the TKE budget across the wall-normal direction, though not in a monotonic and consistent way. The mechanisms influencing the production and the dissipation of turbulent kinetic energy will be discussed further in the next chapter.

Turbulence coherent structures
It is commonly accepted in the wall turbulence research community that turbulence is organized into structures over several length and time scales (Eaton, 1994). According to Adrian (2007), wall turbulence is characterized by the presence of packets of hairpin vortices and their associated quasi-streamwise vortices (QSV) near the wall. The term coherent highlights the fact that these structures possess a "long temporal persistence" in the flow, and their existence has been extensively studied both experimentally and numerically. In particular, a mechanism to generate QSV in a turbulent channel flow has been suggested by Zhou et al. (1999). In their paper, the authors argue that the turbulent boundary layer contains a large number of hairpin vortices, aligned in the streamwise direction as coherent packets. Though they first studied the evolution of a single ideal symmetric hairpin, they went on to investigate asymmetric hairpins, and concluded that the QSV generated from asymmetric structures occurred most often singly and rarely as counter-rotating pairs of equal strength. Zhou et al. added that it was more likely to observe the asymmetric vortices generation process in natural turbulent boundary layers observed experimentally. On the other hand, Jeong et al. (1997) investigated the role of the coherent structures near the wall in a turbulent channel flow, and presented a model of overlapping and alternating-sign QSV as the dominant near-wall structure. Figure 15: Schema of an array of QSV, side view (Jeong et al., 1997) 58 The figure 15 illustrates schematically the disposition of the QSV at the wall.
Again, the convention used for the coordinates is such that the y-direction corresponds to the wall-normal direction here, in a similar fashion as KMM (1987).
The The left-hand side corresponds to the case without any particles, while the right-hand side shows u and ω y for the case IV, representative of the several particle-laden turbulent flows. In our study, the vorticity was computed via a least-square numerical scheme, from the general expression of the vorticity: , with ijk being the alternating tensor.
In particular for channel turbulence where the flow is strongly anisotropic, the vorticity component of interest is: The numerical expression then becomes: computed at a given time, once the statistically steady state has been reached.
It is clear that these snapshots capture the dominant turbulence structure occurring in our computational domain. In particular, the pattern of vorticity Conversely, the velocity is maximal near the center of the channel.
The upper right plot exhibits the striking absence of structures at the walls, as well as a reduction of the velocity at the centerline of the channel. The vorticity ω y is also significantly receded, except a few traces of negative vorticity at the walls. It is evident that the major effect of small inertial particles is to dampen the turbulence coherent structures near the boundaries of the channel. The lower row of instantaneous snapshots indicates that the effect of the particles on the turbulence structure is less obvious, other than the fact that the mean velocity is slightly reduced.

Summary
First, we have validated the implementation of the turbulent GPGPU code to simulate turbulent boundary layer flow without particles against the reference paper of Kim et al. (1987). We have concluded that despite a smaller computa- In the clean channel case, the Reynolds and the viscous forces compensate each other exactly, so that their sum equals the mean pressure gradient applied to sustain the Poiseuille flow. For the cases I to VI, the particle feedback appears in the force balance, and a new force distribution 62 takes place: both the Reynolds force and the viscous force are much lower than in the clean channel, while the particle force increases progressively for larger mass fractions and decreases for larger Stokes number (since a lower number of particles is required to reach the given mass fraction of 0.25 for cases V and VI) . Beyond the momentum conservation, the turbulent kinetic energy (TKE) budget has been investigated, and the mean-shear production, the dissipation rate and the particle production have been evaluated. The particle production is usually small relative to the production term, and can change sign depending on the particle size (cases I and III: P p < 0, case V: P p > 0, case VI: P p ≈ 0). Both P and ε are found smaller than the baseline in all six cases plotted, implying that the general effect of particles on turbulence is to reduce the production and attenuate the dissipation. Finally, we have investigated the interactions between the turbulence coherent structures and the dispersed phase. From the various snapshots of instantaneous velocity and vorticity, we may argue that the coherent structures (quasi-streamwise vortices and sweeps) are weakened near the channel walls in presence of particles.

Discussion
In this chapter, we first discuss our results in light of the current knowledge of turbulent particle-laden flows, focusing on the coherent structures and the quadrant analysis of the turbulence Reynolds shear stress. Next, we interpret the results within the context of Oceanography and discuss the role of sea spray in modifying the momentum exchange between the atmosphere and the ocean.

Coherent structures-particles interactions
In previous studies of particle-laden flows one notable finding has been a strongly heterogeneous spatial distribution of the particles. This phenomenon has been reported both in numerical studies and laboratory experiments, to various degrees. This behavior does not seem to depend on the coupling regime considered: the tendency of particles to segregate and accumulate near the wall remains consistent whether the particle-laden flow is in the one-way, two-way or four-way coupling regime.
In order to discuss this phenomenon in more detail, let us first introduce the quadrant analysis of the Reynolds shear stress. The quadrant analysis divides the Reynolds shear stress −ρ f u w into four classes of events according to the signs of u and w (e.g., Willmarth and Lu (1972) and Wallace et al. (1972)): • the first quadrant, with u > 0 and w > 0, • the second quadrant, with u < 0 and w > 0, • the third quadrant, with u < 0 and w < 0, and • the fourth quadrant, with u > 0 and w < 0.
It is common to denote the various quadrants events by Q1, Q2, Q3 and Q4 respectively. In particular, Q1 events are outward motions of high-speed fluid, Q2 events are called ejections, and represent ejections of low-speed fluid away from the wall, Q3 events are inward motions of low-speed fluid, while the last quadrant part (Q4) contains inrushes of high-speed fluid, also named sweeps (KMM, 1987).
The QSV described by Jeong et al. (1997) and Zhou et al. (1999) are ejection events belonging to the second quadrant, with negative streamwise velocity and positive vertical velocity. It is largely accepted by the wall-turbulence community that these quasi-streamwise vortices are supposedly responsible for the existence of low-speed streaks close to the walls. Marchioli and Soldati (2002) investigated the mechanisms of particle transfer and segregation in a turbulent boundary layer through DNS experiments. They defined turbophoresis as the phenomenon of particle migration toward the wall.
The authors suggested that both sweeps and ejections (Q4 and Q2 events) were efficient transfer mechanisms for the dispersed phase, and proposed a scenario to explain the particle transfer in the channel flow. The figure 18 is from their original paper and illustrates the action of the coherent structures on the small particles.
According to this study, particles are transferred by sweeps in the wall region, where they preferentially accumulate in the low-speed streaks environments, whereas ejections transfer particles from the wall to the outer flow. Particles tend to accumulate in the viscous sublayer since the sweeps dominate ejections. Adrian (2007) detailed the organization of coherent structures in wall turbulence, in the absence of inertial particles. He explained how the low-speed regions observed near the wall were associated with QSV lifting fluid upwards away from the wall. Besides the study of Marchioli and Soldati (2002), there have Figure 18: Schema of particle transfer, From Soldati and Marchioli (2009) been numerous papers studying the coupling between small particles and wall turbulence structures. Nowadays, it is commonly accepted that inertial particles are mainly found in low speed, low vorticity and high strain rate regions, near the channel walls, and that the QSV tend to transfer the particles from the wall to the outer flow (Li et al., 2012).
Another aspect of the turbulence modulation is related to the inertia that the particles display. Since particles have a finite inertia (or a finite particle response time τ p and Stokes number St), they do not exactly follow the turbulent eddies. As a result, they tend to break the eddies and reduce the spatial scale of the coherent structures in the channel -sometimes to the point of eliminating, as observed in figure 16. In complement to the instantaneous plot of the vorticity ω y , we have 66 computed the probability of the four quadrant events for the clean channel and a representative particle-laden flow, at two different heights and have plotted it in figure 19 (see Li et al., 2012). Not surprisingly, the probability of ejections and sweeps at the wall in the clean channel are the largest, confirming that they are indeed the dominant coherent events at the boundary. At the mid-height of the channel, the distribution of the various Q events becomes more uniform, though there are more Q4 events than the other three types of quadrant events. Another important particle effect appears in the reduction of the production term P in the TKE budget. Recall that P is proportional to the product of the Reynolds shear stress and the mean shear. KMM (1987) described the second-and fourth quadrant events as contributing to positive mean-shear production in the budget of TKE. By breaking down the vortical structures, the particles reduce effectively the Reynolds stress (as seen in figure 13). This then leads to the reduction of P . The decrease of the wall-normal component of the r.m.s. of the velocity for all cases (to various degrees) is also consistent with the reduced turbulent stress.

Sea-spray problem
In this study the sea spray droplets present at the sea surface were represented by small solid and inertial particles dispersed in a turbulent Poiseuille flow. Within the framework of this project, we were interested in the effects of the size and concentration on the flow. The figure 20 summarizes the several DNS runs we performed, in comparison to the studies of  and Zaho et al. (2013). On the parameter space St + −φ, the cases I to VII have been plotted as small red crosses, while the data points from Zaho et al. (2013) are represented by small blue circles and  by small blue squares. Our strategy was to cover parts of the two-way coupling region by varying progressively φ (cases II, III and IV), and then increasing St + (cases I, III, V, VI). The case VII was performed to observe the possible impact of very large particles at high concentration. However, the most of the results concerning the cases VII were not shown here, since they were not significantly different from the results of case VI.
In their latest paper,  mentioned that the Kolmogorov length scale η K (defined in chapter 2) was of the order of 1mm in the case of real sea surface during high wind events. It is mostly accepted in the Oceanography community that sea spray is a polydispersed system composed of droplets of various sizes. The size spectrum ranges from tens of microns for the smallest droplets to a few millimeters for the largest spume droplets (Veron, 2012).
Accordingly, the range of Kolmogorov-based Stokes number values goes from 0.2 The determination of the concentration of sea spray has been proven to be challenging: in-situ measurements are scarce, and laboratory experiments are constrained by the size of waves (up to a few meters for the waves generating the droplets). It is likely that the actual sea spray concentration varies over a large range depending on wind speed and sea states. We have therefore focused on the concentration range of the two-way coupling in this study.
We retrieved the Kolmogorov Stokes number for all seven particle-laden cases in our numerical studies and plotted its values in the figure 21. As discussed previously, St K varies across the channel height due to the variation of the dissipation rate ε from the wall to the center of the channel where it becomes very small (see figure 14). The Stokes number St K varies from around 220 for cases VI and VII While the original sea spray problem was greatly idealized, some conclusions can be drawn from the DNS experiments. The main role of inertial particles is to significantly modify the turbulence structure near the wall. Specifically, particles tend to accumulate near the wall and break down the turbulent coherent structures present in the turbulent boundary layer, via the cross-trajectory effect.
Consequently, the Reynolds stress and the mean-shear production of turbulence are reduced, and the spanwise and the wall-normal turbulence intensities are also reduced. The streamwise root-mean-square fluctuations are enhanced, however, due to the presence of particles, which agrees with results found in other studies (Eaton et al., 1994).
The impact of the particles on the mean velocity profile is less clear. While the turbulence structure is significantly modified near the wall, the mean velocity profile is hardly affected near the wall. Although the mean velocity in the middle of the channel is consistently reduced in our numerical experiments, the impact is quite small. Furthermore, Zaho et al. (2013) report an opposite effect (increase of the mean velocity) in their numerical simulations. Clearly, further detailed studies are needed to investigate the effect of particles on the mean velocity profile.

Conclusions
The aim of this project has been to investigate the impact of sea spray droplets on the air turbulence right above the sea surface and the air-sea momentum transfer at strong winds. As a first step to study this complex problem, we have investigated a turbulent particle-laden flow, with the strong winds playing the role of the carrier fluid (gas phase), and the small sea spray droplets being represented by pointwise inertial and solid particles of the same density (dispersed phase).
By adopting the Eulerian-Lagrangian approach to model the coupling between the fluid phase and the particles, we have saved considerable computational time and resources and have simulated a turbulent Poiseuille flow where particles are dispersed and allowed to exchange momentum with the carrier phase.
The mathematical framework of our project consists of the particle equation of motion and the Navier-Stokes equations, plus a feedback term to take into account the effects of the particles on the fluid (Elghobashi, 2003). The equation of motion is applied here under its most reduced form, after considering the magnitude of the various forces present in the complete equation (Maxey and Riley, 1982), discussed in chapter 2. The fluid governing equations are solved through the innovative Lattice Boltzmann method, in complement to the particle equation.
Despite being limited in our computational capacity with only one GPGPU available, we may reproduce the turbulent Poiseuille flow that is in a good agreement with the benchmark study of Kim et al. (1987). The mean velocity, the Reynolds stresses, as well as the two-point correlations, the skewness and the flatness of the velocity fluctuations are all close to the original results. After validating the implementation of the turbulent flow simulation, the particles have been introduced to the channel flow. Using the results we have evaluated the effect of particles on various quantities, such as the Reynolds and viscous forces, the production or the dissipation rate for the turbulent kinetic energy. The computational limitations we faced during this project can be dealt with by using a multi-GPUs approach, or even a non-uniform grid with a finer resolution at the walls and a coarser resolution at the center of the channel. Indeed, the DNS of a turbulent channel flow using two GPUs shows that the statistics of the non-laden case converge better to the results obtained by KMM, since the computational domain becomes larger (Banari, personal communication, 2013). On the other hand, using a single GPU with a non-uniform grid prevents from overresolving the flow in the horizontal directions while getting a higher resolution at the wall, hence saving computational memory.
Our numerical results are mostly in good agreement with the findings of the previous studies on the turbulent particle-laden flows. The streamwise fluctuations (r.m.s.) are increased, whereas the spanwise and the wall-normal turbulence intensities are damped by the action of the particles. As a result, the Reynolds force and the viscous force are reduced in the force balance. The mean-shear production and the dissipation are found smaller than in the case of the clean wall turbulence. The break down of the vortical structures at the boundary by the particles is a likely reason why the Reynolds stress and the shear production of the turbulent kinetic energy are diminished.
In summary, our study suggests that the sea spray droplets suspended in an airflow just above the air-sea interface can significantly modify the near surface 73 turbulence characteristics. However, the effect of the sprays on the drag coefficient (air-sea momentum flux) remains unclear. Our results show slight decrease of mean streamwise velocities at the channel center due to the particles, in contrast with the results of Zaho et al. (2013), which showed slight increase of the mean velocities.