The Impacts of Gustiness on Air–Sea Momentum Flux

: The exchange of momentum across the air–sea boundary is an integral component of the earth system and its parametrization is essential for climate and weather models. This study focuses on the impact of gustiness on the momentum ﬂux using three months of direct ﬂux observations from a moored surface buoy. Gustiness, which quantiﬁes the ﬂuctuations of wind speed and direction, is shown to impact air–sea momentum ﬂuxes. First, we put forward a new gustiness formula that simultaneously evaluates the impact of ﬂuctuations in wind direction and speed. A critical threshold is established using a cumulative density function to classify runs as either gusty or non-gusty. We ﬁnd that, during runs classiﬁed as gusty, the aerodynamic drag coefﬁcient is increased up to 57% when compared to their non-gusty counterparts. This is caused by a correlated increase in vertical ﬂuctuations during gusty conditions and explains variability in the drag coefﬁcient for wind speeds up to 20 m/s. This increase in energy is connected with horizontal ﬂuctuations through turbulent interactions between peaks in the turbulent spectra coincident with peaks in the wave spectra. We discus two potential mechanistic explanations. The results of this study will help improve the representation of gustiness in momentum ﬂux parameterizations leading to more accurate ocean models.


Introduction
In the marine boundary layer, the ocean is tightly coupled to the atmosphere through the wind stress τ. The wind stress is a quantification of the momentum flux across the air-sea interface which is important for many processes such as wave growth and breaking [1][2][3], aerosol production [4], atmospheric and oceanic circulation (e.g., [5][6][7]), global climate [8], and upper ocean mixing (e.g., [9,10]). The momentum flux is also a key component of tropical cyclone intensity [11].
In the constant-flux layer the total wind stress can be decomposed into three parts: Here, τ t is the turbulent stress, τ w the wave-induced stress, and τ v the viscous stress [12]. τ v is caused by the differential of wind and sea surface current speed within the viscous sub-layer which has a thickness on the order of 10 −5 m, much thinner than the wave-induced boundary layer [12,13], and is negligible away from the sea surface. τ w occurs in the wave-boundary layer, within which air flow is influenced by waves [14]. The wave-boundary layer height has been shown to depend on wind speed [15,16] and sea state [17], and is predicted to be order 1 m for pure wind sea conditions in [18]. It has also been demonstrated that the wave induced stress is less than 5% of the total stress at the height of 1 m, and less than 1% at 3 m. Previous experiments have established the practice of making measurements 3-4 m above the sea surface e.g., [19,20] where viscous and wave-induced stress were deemed negligible. (i.e., τ ≈ τ t τ w τ v ). Under such conditions, turbulent fluxes approximate the total momentum flux, e.g., [21]: Here u , v , and w , are velocity in the along, cross, and downwind directions, respectively, with primes denoting the fluctuating component and overbar representing time average (O~30 min). Air-sea coupling within the marine boundary layer is widely parameterized by the wind stress equation: where u * is the shear or friction velocity of the wind at the sea surface, ρ is the air density, C d is the drag coefficient, and U 10 is the wind speed at a reference height of 10 m. Combining Equations (2) and (3), we arrive at C d , calculated from direct measurements Because of the isolation and harsh environmental conditions, measuring fluxes at sea is very challenging and expensive, especially at high wind speeds. As such, it is common for the momentum flux to be determined instead by using the nondimensional drag coefficient as a function of the mean wind speed. Since the drag coefficient is widely adopted and used in lieu of direct flux measurements when quantifying air-sea momentum exchange (e.g., [22,23]), we focus on the impact of gustiness on the drag coefficient such that results have direct implications for Earth system models.
Many studies have quantified the influence of sea state on C d . Results show that wave age C p /U 10 , where C p is the phase speed of the peak wave frequency, is well correlated to C d [24,27,32,33], and C d tends to be higher in young seas, C p /U 10 < 1.5. The association between wave age and C d is not surprising given that young waves are steeper and rougher than old waves [29,34,35]. Also based on the association between wave age and drag, Oost et al. [36] found a reliable relationship between for C p and ocean roughness.
The impacts of swell on C d have also been extensively investigated with mixed results. For example, Dobson et al. [37] found no evidence that swell influences C d , but Drennan et al. [38] observed reduced wind stress in the presence of swell. This finding was supported by the finds of Potter [39], who determined that swell can reduce C d up to 37% due to the reduced turbulent energy around the swell frequency. In contrast, Drennan et al. [40] found that both crossing-and opposing-swell significantly increase C d , and Vincent et al. [41] showed an increase in drag due to increased energy in the high-frequency waves relative to pure wind seas, consistent with their previous work [42]. Increased drag coefficient in the left-rear of a tropical cyclone where conditions are thought to be predominantly crossing and opposing-swell, has also been shown [28].
Surface currents have been shown to reduce C d , especially when the current is aligned with the wind [30]. Kara, Metzger, and Bourassa [7] reported that C d can be reduced over 10% in the tropical Pacific Ocean due to currents being in line with the prevailing winds. Furthermore, a numerical modelling study of an idealized typhon showed that surface current can reduce the momentum fluxes in the right-rear typhoon quadrant about 6% in [5]. In contrast, Wüest and Lorke [43] found that currents in the Kuroshio Extension have very little influence on C d . Currents can also affect the momentum fluxes in indirectly; Kenyon and Sheres [44] found that strong current gradients impact the surface gravity wave field which can then affect air-sea momentum fluxes. The stability of atmosphere influences C d according to Monin-Obukhov similarity theory [45][46][47][48][49]. Monin and Obukhov [50] put forward the Obukhov length to describe the air-sea fluxes processes.
where k is the von Karman constant, T 0 is surface temperature, g is gravity acceleration, v * is friction velocity, q is kinematic heat flux, C t is specific heat, and ρ is air density. The Obukhov length accounts for the relationship between sea surface temperature, specific heat, air density, heat flux, and friction velocity. In light winds, the momentum fluxes are sensitive to the stability of the atmosphere, with increased flux for unstable conditions [51]. While much attention has been given to sea state, currents, and atmospheric stability, there are few studies that directly address the effect of gustiness on C d . The importance of accounting for gustiness when estimating momentum fluxes was first addressed by Janssen [52], who believed that the friction velocity should be related to both the mean and standard deviation of wind speed. Miles and Ierley [53], using a theoretical model, found that, under some circumstances, 20-30% more energy is transferred to the ocean during gusty conditions but that a reduction in the stress is also possible. Abdalla and Cavaleri [54] used a model to compare the wave field under normal and gusty winds and found that wave height is increased under gusty conditions, implying increased momentum flux. Ponce and Ocampo-Torres [55] also showed that wind speed variability increases wave height while Pleskachevsky et al. [56] remotely observed higher wave height under the influences of gustiness. Annenkov and Shrira [57] studied wave responses to different gustiness regimes using an idealized numerical model. They found that gustiness can impact growth rate during wave development but had little impact as waves approached full development. Uz et al. [58] found that decreasing winds had higher wind stress than the increasing winds for any given wind speed, which is mainly caused by the delayed response of short wind waves to the modulated wind forcing.
Previous research notwithstanding, it is very difficult to account for the impact of gustiness on momentum flux because typical wind speed data, whether extracted from satellite remote sensing, collected in situ, or model derived, are produced as an average over space and/or time and so high frequency fluctuations are unknown. As such, it is difficult to incorporate gustiness information into C d parameterizations. Babanin and Makin [31] addressed this by developing a method to estimate maximum gustiness using U 10 , which can be then used to estimate the contribution of gustiness to C d for any given wind speed. In their study, gustiness was quantified using the mean and standard deviation of wind speed calculated over 30 min: Following Babanin and Makin [59], Ting et al. [60] found that gustiness can increase scatter in C d , and used a scale analysis based on the equation of motion to determine that wind stress was proportional to gustiness. They speculated that the variability in the drag coefficient was caused by gustiness through distortion of the air-sea boundary layer structure. The Coupled Ocean-Atmosphere Response Experiment (COARE) algorithm [30] represents gustiness as boundary layer-scale eddies and is more akin to vertical thermal stability rather than horizontal wind fluctuations.
In this study, we focus on the influence of gustiness on the momentum flux. We introduce a gustiness formula that is sensitive to the fluctuations of both wind speed and direction. The metric captures variation within scales of 0.01-1 Hz. We apply the new formula to data from the Impact of Typhoons on the Ocean in the Pacific (ITOP) experiment [61]. Our results indicate that gusty periods exhibit higher C d than predicted by the linear wind speed parameterization [25]. We also used a scaled Miyake-Cospectrum [62] to explore the momentum fluxes transferred into the ocean. This Miyake-Cospectrum is normalized by wind speed, which is helpful for comparing the gustiness influences under different mean wind speeds. We ultimately conclude that gustiness is important to C d , accounting for about 57% of the variability in C d . We also find that the cospectra of the momentum flux is proportional to the gustiness. When winds are gusty, we find an associated increase in energy in the cospectra.

Data
Data used in this study were collected during the Impact of Typhoons on the Ocean in the Pacific (ITOP) experiment [61]. Data were collected between August and December 2010, from an instrumented Extreme Air-Sea Interaction (EASI) buoy (Drennan et al. [63]) deployed in the Philippine Sea approximately 750 km east of Taiwan. Direct fluxes were measured onboard with simultaneous mean parameters including air and ocean temperature, humidity, and atmospheric pressure. The EASI buoy, which was based on the hull of a 6-m Navy Oceanographic Meteorological Automatic Device (NOMAD), also measured directional wave spectra that were produced over 30 min, as was wind speed [64]. Several studies have previously used data collected by EASI during ITOP (e.g., [9,27,39,63,65]) and provide extensive information about the experiment, instrumentation, and data processing and quality. Here, we only reproduce information about the wind sensors and associated processing methods.
The EASI buoy was equipped with two Gill R2A sonic anemometers. However, due to severe weather, only one R2A anemometer collected reliable data. The sampling frequency of the wind anemometer was 20 Hz with runs of 30 min, which was found suitable for collecting all turbulent scales [66]. The anemometer was positioned 5.45 m above mean sea level and away from the buoy structure such that flow distortion was negligible. Two motion packages in the buoy measured six degrees of freedom of motion along with the compass bearing. The motion of the platform was accounted for when calculating the wind velocity in accordance with the methods by [67]. Following motion-correction, a quality control procedure removed spikes in the 20 Hz u, v, and w time series [66]. Ogives were calculated and it was revealed that the scales of 0.01-1.0 Hz were adequate to fully characterize momentum flux, as expected from universal curves [62]. Next, u, v, and w were rotated according to the mean wind vector so that the w was around 0 and u pointed towards the mean wind direction, thereby mean(v) = 0. These data were then used to calculate C d following Equation 4. Log layer theory was used to adjust the value of C d to the equivalent measured at 10 m in neutral conditions, which is referred to as C d10N .

Gustiness Formulation
The Babanin and Makin [31] equation for gustiness (Equation (5)) quantifies the fluctuation of wind speed but does not consider wind direction variability. We show that directional variability can also increase variance in both u and v, and therefore must be part of the gustiness formulation in order to more fully characterize wind variability. For example, if wind speed is constant, the runs with unstable direction will contain more turbulent energy than runs with stable direction. As such, we put forward a new gustiness formula, G A , (Equation (7)), to account for fluctuations in both wind speed and direction.
Equation (7) can be expanded as below: The cov(u , v ) term represents the influence of changing wind direction, and it reduces to 0 when the wind direction is stable. Supposing at one moment that the mean wind direction has turned about x degrees, it will add an extra fluctuation (u , v ) ∼ U·(1 − cosx, sinx) to the original wind vector. A change in wind direction projects an extra velocity onto the (u , v ). It can be seen that the added fluctuation (U − Ucosx, Usinx) is correlated, which is one source of the covariance term. One thing to note, cov(u , v ) can be underestimated to some extent given that some portions of (u , v ) are still independent, similar to Brownian motion [68], when the wind direction is changing. When cov(u , v ) = 0, observation samples of (u , v ) are distributed symmetrically but the main axis is not along the x or y axis. Graphical representations of cov(u , v ) = 0 and cov(u , v ) = 0 are shown in Figure 1. Equation (7) can be expanded as below: The ( , ′) term represents the influence of changing wind direction, and it reduces to 0 when the wind direction is stable. Supposing at one moment that the mean wind direction has turned about x degrees, it will add an extra fluctuation ( ', ')~ • (1 − , ) to the original wind vector. A change in wind direction projects an extra velocity onto the ( , ). It can be seen that the added fluctuation ( −  ,  ) is correlated, which is one source of the covariance term. One thing to note, ( , ) can be underestimated to some extent given that some portions of ( , ) are still independent, similar to Brownian motion [68], when the wind direction is changing. When ( , ) ≠ 0, observation samples of ( , ) are distributed symmetrically but the main axis is not along the x or y axis. Graphical representations of ( , ) = 0 and ( , ) ≠ 0 are shown in Figure 1. To better characterize , sensitivity to sampling frequency was explored. Figure 2 shows the standard deviation of ' + ' for an arbitrary run as a function of sampling frequency. When the frequency reaches approximately 1 Hz, ( ' + ′) converges upon a constant value, in this case 1.83. This is a typical result (though other runs converge upon different values), indicating that our sampling frequency of 20 Hz is more than adequate to obtain a stable standard deviation of and ′.
Our definition of gustiness is sensitive to run length because of the increased potential for wind direction change with time. Thirty-minute run-time is a compromise between a length adequate to capture all turbulent length scales (e.g., [19,37,63]) and short enough such that stationarity is not violated. This record length was used in many previous studies focused on Cd (e.g., [19,21,24,38,66,69]), and, furthermore, 30 min is a typical time domain used by air-sea coupled climate or hurricane models (e.g., [70,71]). Hence, by To better characterize G A , sensitivity to sampling frequency was explored. Figure 2 shows the standard deviation of u + v for an arbitrary run as a function of sampling frequency. When the frequency reaches approximately 1 Hz, std(u + v ) converges upon a constant value, in this case 1.83. This is a typical result (though other runs converge upon different values), indicating that our sampling frequency of 20 Hz is more than adequate to obtain a stable standard deviation of u and v .   Figure 3 shows the times series of wind speed ( ) and gustiness . High can be attributed to the fluctuation of wind speed and direction or low , which increases in Equation (7). Our definition of gustiness is sensitive to run length because of the increased potential for wind direction change with time. Thirty-minute run-time is a compromise between a length adequate to capture all turbulent length scales (e.g., [19,37,63]) and short enough such that stationarity is not violated. This record length was used in many previous studies focused on C d (e.g., [19,21,24,38,66,69]), and, furthermore, 30 min is a typical time domain used by air-sea coupled climate or hurricane models (e.g., [70,71]). Hence, by adopting a 30 min sampling period, we seek to maximize the applicability of the results to future research. Figure 3 shows the times series of wind speed (U 10N ) and gustiness G A . High G A can be attributed to the fluctuation of wind speed and direction or low U 10N , which increases G A in Equation (7).  Figure 3 shows the times series of wind speed ( ) and gustiness . High can be attributed to the fluctuation of wind speed and direction or low , which increases in Equation (7).  Figure 4 shows the cumulative probability function of for all the runs to help define a limit above which the wind is gusty.

Results
of gusty and not-gusty runs is calculated via Equation (6). A value of 0.165 was chosen for selecting gusty runs and represents the time when the CPF is 50%.  Figure 4 shows the cumulative probability function of G A for all the runs to help define a limit above which the wind is gusty. G A of gusty and not-gusty runs is calculated via Equation (6). A value of 0.165 was chosen for selecting gusty runs and represents the time when the CPF is 50%.  Based on the formula for , Figure 5 shows as a function of . Each data point is color coded by gustiness, with dark green indicating 0.165. When the winds are not gusty, a clear linear relationship between & emerges. This relationship is similar to the linear empirical equation proposed by [25,30]. This suggests that these widely used linear equations were established based on data collected during relatively wind speed stable conditions. Most gusty points are over this linear relationship, suggesting that gusty conditions are generally associated with an increase in momentum flux. From the error bar plot, for moderate winds (5-12 m/s), we find that the variation in of winds can be reduced by as much as 57% by removing the gusty runs. Based on the formula for G A , Figure 5 shows C d10N as a function of U 10N . Each data point is color coded by gustiness, with dark green indicating G A ≤ 0.165. When the winds are not gusty, a clear linear relationship between C d10N & U 10N emerges. This relationship is similar to the linear empirical equation proposed by [25,30]. This suggests that these widely used linear equations were established based on data collected during relatively wind speed stable conditions. Most gusty points are over this linear relationship, suggesting that gusty conditions are generally associated with an increase in momentum flux. From the error bar plot, for moderate winds (5-12 m/s), we find that the variation in C d10N of winds can be reduced by as much as 57% by removing the gusty runs. emerges. This relationship is similar to the linear empirical equation proposed by [25,30]. This suggests that these widely used linear equations were established based on data collected during relatively wind speed stable conditions. Most gusty points are over this linear relationship, suggesting that gusty conditions are generally associated with an increase in momentum flux. From the error bar plot, for moderate winds (5-12 m/s), we find that the variation in of winds can be reduced by as much as 57% by removing the gusty runs.  The covariance term will be increased when the direction is turning. Therefore, Figure 6 indicates that changing wind direction could play an important role in contributing to the gustiness, and thus . In another words, some scatter in , or increased gustiness, can be explained by changes in direction. To parallel this viewpoint, the data are plotted again (Figure 7), but are now identified based on the stability Figure 5. C d of gusty winds (red & yellow points) and not-gusty winds (green points) selected by G A . Non-gusty winds are selected by the condition G A ≤ 0.165 and marked as green. The blue line is the averaged C d for G A ≤ 0.165 with error bars denoting 95% confidence interval binned every 2.5 m/s. The orange line is the averaged C d for G A > 0.165 with error bars also denoting 95% confidence interval binned every 2.5 m/s. The black dash line is the C d10N predicted by [25]. The red line is C d10N predicted by COARE 3.0 algorithm [30]. Figure 6 shows C d10N vs. U 10N . Points are colored by the covariance term, |cov(u , v )|/U 10 . This term can be attributed to the changing wind direction, since a change in direction (θ) will add extra correlated oscillation u ∼ U 10 ·(1 − cosθ), v ∼ U 10 ·sinθ to the fluctuation. The covariance term will be increased when the direction is turning. Therefore, Figure 6 indicates that changing wind direction could play an important role in contributing to the gustiness, and thus C d . In another words, some scatter in C d , or increased gustiness, can be explained by changes in direction. To parallel this viewpoint, the data are plotted again (Figure 7), but are now identified based on the stability of their wind speed and direction. Here, we used std(u) > 1m/s and std(θ) > 7 • to identify runs that were considered either wind speed unstable or wind direction unstable. Using these values results in 25% and 31% of data characterized as gusty due to variability in speed and direction, respectively. For comparison, if we set the threshold to 0.5, 2, or 3 m/s, the percentage of gusty data would be 78%, 4%, 0.7%, respectively. Setting the threshold to 3, 5, or 10 • would result in 100%, 82%, and 11%, respectively, classed as gusty. Most importantly, Figure 7 shows that deviations in both speed and direction are important considerations when quantifying gustiness. Notably, directional deviations identify much of the scatter in the U 10N − C d10N relationship at low to moderate wind speeds that cannot be captured when considering fluctuations in wind speed alone. Furthermore, we avoid another issue with identifying gustiness based on wind speed, which is that all data above~15 m s −1 are classified as gusty because, inherently, std(u) increases with U. of their wind speed and direction. Here, we used ( ) 1m/s and ( ) 7° to identify runs that were considered either wind speed unstable or wind direction unstable. Using these values results in 25% and 31% of data characterized as gusty due to variability in speed and direction, respectively. For comparison, if we set the threshold to 0.5, 2, or 3 m/s, the percentage of gusty data would be 78%, 4%, 0.7%, respectively. Setting the threshold to 3, 5, or 10° would result in 100%, 82%, and 11%, respectively, classed as gusty. Most importantly, Figure 7 shows that deviations in both speed and direction are important considerations when quantifying gustiness. Notably, directional deviations identify much of the scatter in the − relationship at low to moderate wind speeds that cannot be captured when considering fluctuations in wind speed alone. Furthermore, we avoid another issue with identifying gustiness based on wind speed, which is that all data above ~15 m s −1 are classified as gusty because, inherently, ( ) increases with .     We find that variations in both wind speed and direction can have significant influences on . In fact, the stability of wind direction is potentially a better metric than the stability of wind speed (Figure 7). When wind direction is stable, a linear relationship suggested by previous studies (e.g., [21,24,66]) between and emerges. In contrast to ( ′, ′), speed stability does not correlate with variation in at moderate wind speeds (5-12 m/s). This implies that direction is more effective in explaining the variation in . We find that variations in both wind speed and direction can have significant influences on C d . In fact, the stability of wind direction is potentially a better metric than the stability of wind speed (Figure 7). When wind direction is stable, a linear relationship suggested by previous studies (e.g., [21,24,66]) between C d10N and U 10N emerges. In contrast to cov(u , v ), speed stability does not correlate with variation in C d10N at moderate wind speeds (5-12 m/s). This implies that direction is more effective in explaining the variation in C d10N .

Gustiness and Vertical Oscillation
Next, we examine in detail the difference between a gusty run and a not-gusty run.  Figure 8a, the data are distributed symmetrically, forming a circle, and the histogram is symmetric, which indicates the wind direction remained relatively stable. In Figure 8b, the data are asymmetrically distributed over a larger area in u , v space, skewed towards the lower-left. In this case, the distribution is an ellipse and the major axis is rotated away from the wind direction, indicating that the direction was changing over the duration of the run. In Figure 8c, the directions of the gusty and non-gusty runs are shown. We can see that the direction of the run in Figure 8b has changed about 20 degrees over 30 min, while the direction of the not-gusty run in Figure 8a has a negligible trend.
Note that the mean of |w | in Figure 8a,b is larger during the gusty run than the nongusty. This suggests a stronger vertical oscillation in the gusty runs which may account for the observed higher C d . To explore this further, we plot the mean of |w | in Figure 9 as a function of wind speed for all our data. This shows that gusty periods tend to have higher |w | compared to periods when the wind speed and direction remain more stable, as determined using G A = 0.165. This indicates a stronger oscillation in the vertical direction during gusty conditions which, when correlated with the u , v , increases C d . For moderate winds speeds, |w | is 30~50% larger for gusty than non-gusty periods. Hence, higher gustiness can inspire a stronger vertical oscillation, which transports the extra energy across the air-sea interface. mained relatively stable. In Figure 8b, the data are asymmetrically distributed over a larger area in ′, ′ space, skewed towards the lower-left. In this case, the distribution is an ellipse and the major axis is rotated away from the wind direction, indicating that the direction was changing over the duration of the run. In Figure 8c, the directions of the gusty and non-gusty runs are shown. We can see that the direction of the run in Figure 8b has changed about 20 degrees over 30 min, while the direction of the not-gusty run in Figure 8a has a negligible trend. Note that the mean of | | in Figure 8a,b is larger during the gusty run than the nongusty. This suggests a stronger vertical oscillation in the gusty runs which may account for the observed higher . To explore this further, we plot the mean of | | in Figure 9 as a function of wind speed for all our data. This shows that gusty periods tend to have higher | | compared to periods when the wind speed and direction remain more stable, as determined using = 0.165. This indicates a stronger oscillation in the vertical direction during gusty conditions which, when correlated with the ', ', increases . For moderate winds speeds, | | is 30~50% larger for gusty than non-gusty periods. Hence, higher gustiness can inspire a stronger vertical oscillation, which transports the extra energy across the air-sea interface. We now offer two potential explanations for the source ( , ) term. The first source is that the wind direction is changing while the other is that the waves may be intermittently coupling with the gustiness. When the wind stress is between wave and wind directions, waves can explain some gustiness. These two sources may function simultaneously. We inspected many cases finding that, for a few runs, the stress direction was in between the wave direction and wind direction.  Figure 10 shows the mean power spectra of ′ for the entire data set. Each subplot represents a different wind speed range, data are further subdivided based on gusty and non-gusty conditions. Concurrent mean wave spectra are overlaid on the turbulence spectra. Turbulent power spectra can be seen to increase with wind speed, as expected; how- We now offer two potential explanations for the source cov(u , v ) term. The first source is that the wind direction is changing while the other is that the waves may be intermittently coupling with the gustiness. When the wind stress is between wave and wind directions, waves can explain some gustiness. These two sources may function simultaneously. We inspected many cases finding that, for a few runs, the stress direction was in between the wave direction and wind direction. Figure 10 shows the mean power spectra of w for the entire data set. Each subplot represents a different wind speed range, data are further subdivided based on gusty and non-gusty conditions. Concurrent mean wave spectra are overlaid on the turbulence spectra. Turbulent power spectra can be seen to increase with wind speed, as expected; however, in each case, the gusty conditions are always more energetic than the non-gusty runs. For the lowest wind speeds, the gusty case has 20% higher energy, increasing to over 50% at the highest wind speeds. The increases are concentrated around the peak frequency. The peak frequencies of the wave spectra are located around the peaks in the w spectra, which shows some coupled or at least correlated behavior between turbulence and waves. This is true whether or not conditions are gusty. This analysis is repeated for u spectra ( Figure 11) and no difference is found between the u power spectrum of gusty and not-gusty runs (p < 0.05). The same is true of v (not shown), i.e., there is no marked increase in energy in horizontal velocity at ocean wave scales. This suggests that the increased C d during gusty conditions is due to increased energy in the vertical oscillation. In each panel, blue represents the gusty runs and the red represents non-gusty. Note that the limits of the yaxis in each subplot are not the same. Shaded area is 95% confidence interval for the ' spec. In each panel, blue represents the gusty runs and the red represents non-gusty. Note that the limits of the y-axis in each subplot are not the same. Shaded area is 95% confidence interval for the w spec. Figure 10. Mean power spectra of ′ combined with wave spectra for (a) 0 < < 5 m/s, (b) 5 < < 10 m/s, (c) 10 < < 15 m/s, and (d) 15 < < 20 m/s wind speeds. The solid lines are the ' spectra (left y-axis), the dashed lines are wave spectra (right y-axis). In each panel, blue represents the gusty runs and the red represents non-gusty. Note that the limits of the yaxis in each subplot are not the same. Shaded area is 95% confidence interval for the ' spec.  Figure 10 shows that gusty runs coincided with more energetic ' spectra. Wh more energetic ′ is correlated with ', gustiness translates directly to increased mom tum flux. Figure 11 shows the ′ spectra, in which there are less changes in the spec To investigate whether this extra energy is correlated with , we examine a probabi density function (PDF) of normalized momentum flux in Figure 12. The momentum f is normalized by friction velocity (− / * ) where * is estimated using [25] Figure 10 shows that gusty runs coincided with more energetic w spectra. When more energetic w is correlated with u , gustiness translates directly to increased momentum flux. Figure 11 shows the u spectra, in which there are less changes in the spectra. To investigate whether this extra energy is correlated with u , we examine a probability density function (PDF) of normalized momentum flux in Figure 12. The momentum flux is normalized by friction velocity (−uw/u 2 * ) where u 2 * is estimated using [25] and the data are binned by G A . Positive values represent momentum transferred into the ocean and negative values are momentum transferred into the atmosphere. The PDFs have asymmetric shapes skewed towards positive values. This means downward momentum fluxes exceed the upward fluxes. As G A increases, the PDFs flatten and the probability of high and low −u w /u 2 * increases, while the probability of −u w /u 2 * around zero decreases. This shows that more energetic w do indeed increase momentum flux to the ocean.

Possible Reason and Evidence for the Increased Momentum Fluxes
To show how much the net momentum fluxes increased as a result of gustiness, we introduce the u w dimensionless cospectrum, which is scaled using the Miyake formula [62]. This spectrum can contain valuable information about momentum fluxes [30,38,72]. To quantify the net momentum flux, we integrate the u w cospectra along the frequency domain.
Here, S uw is the real part of the cross spectrum of u and w. Figure 13 exhibits the cospectra scaled by mean wind speed U, the measurement height z, and u * . Cospectra are shown for different gusty conditions along with the Miyake universal curve, which is widely used in similar studies [38,39,41]. exceed the upward fluxes. As increases, the PDFs flatten and the probability of high and low − ′ ′/ * increases, while the probability of − ′ ′/ * around zero decreases. This shows that more energetic ' do indeed increase momentum flux to the ocean. To show how much the net momentum fluxes increased as a result of gustiness, we introduce the ′ ′ dimensionless cospectrum, which is scaled using the Miyake formula [62]. This spectrum can contain valuable information about momentum fluxes [30,38,72]. To quantify the net momentum flux, we integrate the ' ' cospectra along the frequency domain.
Here, is the real part of the cross spectrum of and . Figure 13 exhibits the cospectra scaled by mean wind speed , the measurement height , and * . Cospectra are shown for different gusty conditions along with the Miyake universal curve, which is widely used in similar studies [38,39,41]. Figure 13. Observed ′ ′ cospectra normalized using the universal scaling of Miyake, Stewart, and Burling [62] and averaged according to the wind gustiness, shown in multiple colored dash lines. The black dashed line is the Miyake universal curve with values taken from [72]. The red  To show how much the net momentum fluxes increased as a result of gustiness, we introduce the ′ ′ dimensionless cospectrum, which is scaled using the Miyake formula [62]. This spectrum can contain valuable information about momentum fluxes [30,38,72]. To quantify the net momentum flux, we integrate the ' ' cospectra along the frequency domain.
Here, is the real part of the cross spectrum of and . Figure 13 exhibits the cospectra scaled by mean wind speed , the measurement height , and * . Cospectra are shown for different gusty conditions along with the Miyake universal curve, which is widely used in similar studies [38,39,41]. Figure 13. Observed ′ ′ cospectra normalized using the universal scaling of Miyake, Stewart, and Burling [62] and averaged according to the wind gustiness, shown in multiple colored dash lines. The black dashed line is the Miyake universal curve with values taken from [72]. The red Figure 13. Observed u w cospectra S uw normalized using the universal scaling of Miyake, Stewart, and Burling [62] and averaged according to the wind gustiness, shown in multiple colored dash lines. The black dashed line is the Miyake universal curve with values taken from [72]. The red is the averaged wave spectrum of all runs, also normalized using scaling of Miyake, and standardized to a maximum value of 0.3.
In Figure 13, the magnitude of each spectrum increases with increasing gustiness. The cospectra therefore present further evidence that gustier winds have higher momentum fluxes. The energy containing region of the wave spectrum is aligned with that of the cospectra. This overlapping of energetic scales suggests some connection or coupling between the surface waves and the winds. The energy of the vertical turbulence spectra is strengthened in the presence of high gustiness. Although gustiness is characterized by horizontal wind variance, it is clear that when it is gusty there is a coincident increase in the vertical wind variance and at a similar turbulent scale-around 0.10 Hz. The vertical wind variance is connected with the horizontal wind variance due to the turbulent interactions and cascade of energy, which is isotropic. We know that, close to the sea surface, some of the vertical wind variance is coherent with the surface waves, however, direct physical connection remains unclear.
Among the possible explanations, we propose two potential mechanisms for how gustiness can cause the stronger w . First, when the wind direction is rotating, it will cause extra horizontal turbulence and produce additional positive vorticity in the z-direction. Additionally, the ocean surface has a frictional force which can consume this vorticity, which is to say a negative vorticity in the z-direction induced by friction. For example, in a cylinder, extra positive vorticity in the upper layer will be balanced by the negative vorticity held in the bottom layer. In this case, this mechanism causes friction upon the sea surface and may explain the increased w during gusty conditions. Another potential mechanism is that, when the wind speed is increasing, the vertical shear will be strengthened, in turn inspiring a lateral vorticity. For example, when U is increasing, it can create a positive vorticity in the mean wind direction. The only way to consume this extra vorticity is by bottom friction. In this case, small-scale lateral circulation will be produced which can connect the bottom wave surface layer with the upper layer. Ultimately, stronger vertical mixing occurs. However, proof of these two potential mechanisms requires more detailed observations to understand energy pathways and how the extra vorticity is balanced. For example, it is advisable to use the Particle Image Velocity method to take snapshots of the vertical profiles of the winds near the air-sea interface in a laboratory water tank or over the ocean. Measurements made at multiple levels may also help to understand the processes when the gustiness becomes increased.

Summary and Conclusions
A new gustiness formula [Equation (8)] is put forward which quantifies the fluctuations of wind speed and direction. We find that it is important to account for wind direction in explaining variability in the drag coefficient. We find that gustiness can increase C d independent of mean wind speed. The nominal linear relationship between C d and U 10 put forward by previous studies (e.g., [19,24,30], etc.) is found to emerge under non-gusty (or less gusty) conditions. From examining spectra, more energetic vertical turbulence is correlated with gustiness and is perhaps induced by some of the same mechanisms that produce gustiness. The wind cospectra and wave spectra share a peak frequency, suggesting a coupling between the two. From investigation of the cospectra and PDFs, the gusty winds have a clear increase in the normalized momentum fluxes, proving that gustiness results in extra momentum flux under the same wind speed. Overall, gustiness is an important factor for understanding momentum fluxes into the ocean and should be considered in an effective parameterization.