INTRODUCTION Bubbles mediate turbulent fluxes at the air-sea interface. Aerosol fluxes have an intimate dependency on bubbles surfacing and, by bursting, releasing carriers into the atmospheric boundary layer. Whitecaps on the open ocean are the product of turbulence at the sea surface, because of wave-breaking events. As we will see later in the paper, the injected turbulence that produces whitecaps is distributed vertically as well. Whitecaps have the visual aspect of foam patches. For a complete definition, we should relate whitecaps with bubbles by considering whitecaps as clusters of surface bubbles. The bubble size distributions within whitecaps may differ, of course, but a larger size distribution of bubbles is prevalent at the sea surface. Knowledge of whitecapping variability can improve the calculation of aerosol fluxes between ocean and atmosphere , therefore one should consider whitecaps as a primary source for aerosol production. In order to have a proper correlation between aerosol production and whitecaps, one should introduce the concept of whitecap fraction or coverage. The whitecap size area measures the magnitude of turbulence injected by the atmosphere into the ocean. Recent whitecap coverage studies show the roles that wave development, wave-wave interaction, wave-current interaction and wind history (the cause of swells) play in influencing the variations of whitecap coverage in addition to wind speed , , . In this sense, whitecaps strongly relate to energy dissipation of waves, the least known process of wave evolution . There is a direct source for surface foam generation. We are referring at a phenomenon which hasn’t been observed on the open ocean accurate enough. It is described by bubble subsurface dynamics in relation with whitecap dynamics. These bubbles organize in bubble plumes, which feed foam at the sea surface. Bubble plumes variability is crucial in assesing other indirect sources for whitecap generation. In this category are included wind speed at sea surface, wave frequency spectrum and wave amplitude, which can all describe the rate of wave energy dissipation. In this research, we correlate the rate of wave energy dissipation, by considering the parameters just mentioned, with bubble plumes dynamics. Wave energy dissipation can be depicted as a reaction of wave-breaking. It can be quantified at the sea surface, analyzing the size area of whitecap fractions. Wave energy dissipation can be rated at the subsurface level by observing how deep bubbles are injected. Therefore, wave energy dissipation has a tight relation with the range of bubble plumes depth. This study should be an intrusive insight on having a complete picture of the causes that influence the variability of whitecaps, which has a direct effect on aerosol production. In the next sections, some details on why bubbles are important in gas fluxes and aerosol production, existing parametrizations, the work done to realize the observations and data collecting in a high winds regime, and further research ideas are presented. By going right to the direct source and analyzing the dynamics of bubble plumes will help us understand with deep accuracy the mechanisms that drive whitecap variability. It has already been proven that, considering whitecap fraction as function of wind speed only, is not enough for our current understanding. That is why the details in the physics of wave energy dissipation will give us a much clearer view on the sensitiviy in the variations of whitecap fractions. THE IMPORTANTCE OF BUBBLES IN GAS FLUXES AND AEROSOL PRODUCTION Bubbles encapsulate water vapors and other gases. We are interested in analyzing bubbles, because they are a direct cause for transfer velocity of different gases. Aerosols are produced as sea spray particles by breaking the surface tension of individual bubbles. The bubble film layers are fragmentated into small droplets (film droplets). Bubbles usually cary salt and organic particles so they are called “dirty” bubbles. The more “loaded” bubbles are, the more effective they are in generating cloud condensation nucleis. The difficulty is that loaded bubbles should, theoretically, have a higher dragging effect, which comes with a decrease in vertical velocity (upward motion) of the bubbles. Coarser salt or organic particles carried by bubbles have a higher inertia, so the “popping” momentum in launching aerosols into the atmospheric boundary layer is inhibited. Molecular diffusion for subsurface bubbles is an important feature for gas transfer, because it actually reduces the transfer velocity of gases. If the molecular diffusion through bubble walls is higher, then the number of available bubbles at the surface, which can organize in foam patches, is smaller. If the foam patches contain gases with low solubility, they become void fractions. There is a relation between void fractions and interstitial water, which involves the transfer velocity. Interstitial water is the area, at sea surface, unaffected by the disruption of surface tension. Its characteristic si that it surrounds the bubbles, and it has an influence on the total gas transfer rate. For highly insoluble gases, the capacity of the interstitial water may be smaller than that of the bubbles even for small void fractions . This is called restriction of gas transfer in a dense isolated plume “suffocation”. High subsurface bubbles molecular diffusion inhibits the buoyancy of bubbles and, hence, the vertical upward velocity of bubbles. In this sense, the Wolf model (2006) shows that there is a higher gas transfer for low Schmidt number, Sc, therefore subsurface lower molecular diffusion, D, generates more gas available at the surface, hence a higher transfer rate. Schmidt number describes the capacity of molecular difussion through the bubble’s wall. According to Wolf (2006), the Schmidt number is equivalent to molecular difussion D, and it can be found in the following mathematical relation: {K}_{b} = {Sc}^{-x}{\alpha}^{-y}, where Kb is a contribution to the transfer velocity between the main air or gas and water reservoirs and it only depends on the partial equilibrium of a bubble during its lifetime, and x and y are positive constants. Low Schmidt number means higher upward vertical velocity for subsurface bubbles (because of low molecular difussion), therefore more bubbles are available for surfacing. Wolf model assumes two cases for the transfer velocity: one without void fraction and the other with 20 % void fraction. It can be noticed that the transfer velocity is reduced at low Schmidt number with 20 % void fraction present. Correlation between observed gases with different solubilities and bubble size distribution would give an insight on what type of fluxes, from the chemical point of view, are more likely to be enhanced. The presence of void fractions assumes that gas is trapped at the surface, therefore the transfer velocity is suppressed. Considering this phenomenon, one can correlate void fraction with whitecap time decay and, therefore, with transfer velocity in low winds, because in high winds condition, this phenomenon might not be so relevant. The reason is because, during high winds regime, high wind stress fragments the bubbles, so gas transfer is enhanced. PARAMETERIZATIONS AND CORRELATIONS Our interest is to go right to the source of whitecap variability and, therefore, of aerosol fluxes. Other fluxes, such as momentum fluxes especially, have no or little dependency on whitecap variability. One should dynamically couple whitecaps and bubble plumes characteristics in order to explicitly quantify the potential of aerosol production. A range of several parameterizations for whitecap coverage have been performed. Whitecap fraction as a function of wind speed only is the most common, and it takes the form of a power-law relation between W and U₁₀, W = a{U}^b_{10}, with coefficients a and b derived from a set of observations. Because it was clear that W cannot vary with wind speed in a linear way, some autors tried to link W with the cubed wind speed, W = a({U}_{10} + b)^3. In the literature of Monahan and Lu, in 1990, there are two stages in the evolution of whitecaps. In stage A the foam is generated by actively breaking waves, where bubbles are entrained and fragmented inside the breaking wave crest, and is also refered as the acoustically active phase . In stage B whitecaps are subjects to the decaying foam patches, where bubble creation processes cease, therefore the newly formed bubble plume becomes acoustically quiescent and evolves under the influence of turbulent diffusion . During the quiescent phase, one can quantify the decay time of whitecap foam by using a simple exponential model implemented by Callaghan et al. in 2012: A(t) = {A}_{0}\exp\left({\tau} \right), where A(t) is time evolving area of whitecap foam during the decay phase, t is time where t = 0 occurs at the time when the foam patch area is at its maximum value, A₀, and τ is a constant called whitecap foam decay time. According to Philips, in 1985, to obtain the whitecap fraction, one can refer to an average bubble persistence time τbub to obtain W: W = ^{\infty}\tau _{bub}c\Lambda(c)dc, where τbub is the bubble persistence time, Λ(c) is the Philips parameter and c is the forward speed of the beaking wave crest. The Philips parameter is a characterization of breaking wave dynamics in terms of breaking wave speed, and is defined such that Λ(c)dc is the meters of breaking crest per square meter of ocean surface . Another proposal of whitecap parameterization comes from Zhao and Toba, in 2001, namely the non-dimensional Reynolds number RH, which is optimized for wind waves with velocity and length scales represented by u* and wave height Hs, respectively, {R}_{H} = _{*}{H}_{s}}{{\nu}_{w}}, where νw is the kinematic viscosity of water, because the Reynolds number is applied on the ocean. The Reynolds number relates the inertial forces of the fluid to fluid viscosity. If the fluid viscosity is low, the Reynolds number is high, therefore the fluid is more susceptible to turbulence. This shows how much wind stress is needed, for the given kinematic viscosity of water, to generate waves with a critical amplitude at which they can break and produce whitecaps. In this sense, Wolf, in 2005, suggested that whitecap coverage, W, can be related with wind speed and wave height, so that can lead to a direct correlation of whitecape coverage with Reynolds number: W \propto {R}_{H}. Considering the dependency of whitecap fractions on sea state as well, not only on wind speed, is essential in order assess their variability. Swells are are particular sea states that inhibit whitecap production. Because swells move faster than the environmental wind speed, they don’t loose energy through wave energy dissipation (decisive in whitecap production), resulting in attenuated wave-breaking occurence. Swells generate a low level wave-driven supergeostrophic jet that acts to disorganize eddies within turbulent flow at sea surface. This translates into removing friction caused by turbulent stress at sea surface, phenomenon responsible for foam cell rupture within whitecaps , which decreases foam time decay. The amplitude of the swell, according to observations realized by Callaghan in 2012, can be calculated by integrating the wave spectrum between frequencies f₁ and f₂: a = 2_{1}}^{{f}_{2}}S(f)df}, where S(f) is the wave spectral density. Geostrophy relates to a horizontal wind flow, aligned with the isobars (in this case, along the vertical), having no contribution on the vertical component of the wind. Usually, a purely geostrophic horizontal wind at sea surface cannot inject enough turbulence to generate wave breaking. Therefore, a downward ageostrophic component or a shift from the purley geostrophy in the wind at sea surface should be present. The wave-driven supergeostrophic jet subdues the downward ageostrophic component (and enhances the upward ageostrophic component or the upward momentum flux), thus it reduces turbulent injection from atmosphere to ocean. The surface value of this wave-driven jet or stress can be described by integrating the wave-induced stress going into each wave component: \left(0 \right) = ^{\infty}\omega\beta\phi(\omega)d\omega, where ω is the wave angular frequency, ϕ(ω) is the frequency spectrum, and β is the dimensionless wave growth parameter. OBSERVATIONAL MECHANISMS We observed the sea state and fluxes in a high winds regime, during the HIWINGS experiment project. Our region of interest was the North Atlantic ocean, south of Greenland, in the path of strong and frequent cyclones. There was an acute positive NAO (North Atlantic Oscillation) index during that period, namely from the beginning of October to middle November, 2013. Because of this phenomenon, a powerfull jet stream and intense horizontal pressure gradients were present, hence the repetitive cyclones. We succeeded to perform the measurements during winds of over 20 m/s, which makes the data and the future results more valuable, because it is something that probably hasn’t been realized in the past. A notable variety of instruments that measured the sea state and fluxes were used during this field experiment. The scientific crew was partitioned in different teams having specific goals regarding their interest for the data, but a unique goal when it comes to understand the big picture of the air-sea coupled system dynamics. Our specific scientifc goal is to understand the mechanisms of primary source for aerosol production, and hence aerosol turbulent fluxes. We already made the consideration that aerosol turbulent fluxes occur via bubbles, whitecaps, and bubble plumes dynamics. Our main focus is on bubble plume and whitecap physics, but on individual bubble physics as well. Insights about the relation between whitecap coverage variability and bubble plume dynamics will be presented in detail in the next section. Regarding the physics of individual bubbles, our interest is, by using the two acoustical resonators, to obtain the bubble size distribution as a function of depth at a wide range of wind speeds . We will be able to compare the bubble size distribution obtained by using this method with an optical bubble size distribution using the subsurface bubble camera. In order to have a deep understanding of this natural process in a high winds regime, we need to harmonize and correlate the measurements performed using this rich variety of instruments. We used a foam camera to detect whitecap fractions, a submerged bubble camera to explore the optical properties of bubbles within bubble plumes, two resonators to assess the acoustical properties of bubbles within bubble plumes, and a sonar to detect the bubble plumes dynamics. All instruments were set in place along a 11 meters spar buoy. Just under the foam camera, we had wires that could measure the sea state. More detailed informations on the sea state can be generated by the Wave Rider, which was separated from the spar buoy. Other instruments took measurements from the ship’s board, including measurements of aerosol and gas fluxes. RESEARCH APPROACHES AND FURTHER STUDY In this section, our scientific objectives are presented. We remind, as well, the specific instruments which generated the data that we will process and analyze. The first basic objective is to detect the temporal behaviour of bubble plumes. This implies analyzing change of bubble plumes penetration depth in time, as well as bubble plumes duration. In this way, a correlation of bubble plumes duration with whitecap fractions variability and with whitecap time decay, as well, would be an interesting completion at our overall understanding. Thus, we will make a time evolution of bubble plumes per each wave-breaking event, and produce a mean over all observed wave-breaking events. The sonar data will be used to perform these calculations. Next, we will assess how variability in bubble plumes depth affects the whitecap time decay. For instance, we will be able to check if a fast decrease in bubble plumes depth influences the longevity of whitecap fractions. To support this correlation between whitecap fractions time decay and bubble plumes dynamics we will use data collected by the sonar and the foam camera. There was already mentioned that the primary source for whitecap production and variabily is the existence of bubble plumes, therefore our intention is to make direct correlations between bubble plumes and other factors like wind speed and sea state. If we can verify and understand these direct correlations, than we will check the validity of a dynamical relation between bubble plumes and whitecap production. Considering this fact, we want to see how the depth of bubble plumes change with sea state. Thus, we will correlate these changes in bubble plumes depth with different wave states, like wave frequency spectrum and wave amplitude. We want to analyze this bubble plumes dynamics during wind driven waves, which are slower than the environmental wind speed, have higher frequencies and lower amplitudes. For comparison, we will make the analysis during swells as well, which are affected by wind history, are faster than the environmental wind speed, have lower frequencies and higher amplitudes. As a result, we will be able to determine how swells, because of low magnitude wave energy dissipation, shape whitecap coverage variability. Conditions that maximize the effect of swells are weak to moderate wind at sea surface, waves (swells) usually coming from the south , and low frequency waves (around 0.1 Hz). A good way to estimate the effect of swells would be to observe them after the storm passes, during post-storm lower winds. We can perform this by using the wind speed ship data and Wave Rider data. As discussed above, wave energy dissipation is a determinant factor for whitecap production. Going to the direct relation between whitecap coverage and bubble plumes, wave energy dissipation can be considered as a function of bubble plume depth. Significant bubble plume depth translates into large magnitude wave energy dissipation. We must include the whole environmental set when it comes to wave energy dissipation. Therefore, we must consider wind speed, wave frequency and wave height. On the open ocean, winds faster than waves extract energy from the waves by reducing their height, enlarging their slope (having higher frequency), and, thus, leading to wave energy dissipation. The presence of this phenomena sustaines a strong correlation between bubble plumes depth and wave state. Accounting for the fact that wave energy dissipation is proportional to wind speed and bubble plume depth, and inversely proportional to wave amplitude and wave speed (wave phase velocity), we can infer that whitecaps have the same properties. We can show this complex relation of wave energy dissipation with the other parameters through the following equation: $$\epsilon = _{\ast}{d}_{bp}}{{a}_{w} c},$$ where u* is wind speed (wind stress), dbp is bubble plume depth, aw is wave amplitude, and c is wave phase speed. This mathematical relation shows that high wind speed and a deep bubble plume relative to reduced wave amplitude (high wave frequency) and low wave phase speed results in a large wave energy dissipation magnitude. A specific method of estimating the bubble plume depth can be realized by calculating the total air fraction on the vertical. In order to perform this, all bubble sizes, for a given volume, must be integrated. Another interesting and, we think, relevant approach in quantifying whitecap coverage variability is by looking at bubble plumes advections. We want to see if we can grasp differences in the general pattern of bubble plumes advections due to a certain turbulent phenomenon called Stokes drift velocity. This Stokes flow is a turbulent factor within Langmuir circulation and it stands out in the form of horizontal vorticity advection, which is differential, meaning that it changes with depth (in this case, it has a lower magnitude with depth). The Langmuir circulation is the manifestation of relatively shallow turbulence in the upper ocean, and it consists in counter-rotating vortices within a shallow water layer. If the differential vorticity advection due to Stokes flow adds up to general flow pattern of bubble plumes advection, then there is a net shift in the position of bubble plumes relative to the position of whitecaps. Our intention is to observe a sensitivity in whitecap fractions variability due to this shift. This will allow us to find, through the anticorrelation between whitecaps and bubble plumes due to Stokes drift, a direct source of whitecap coverage variability.