The momentum transfer from the atmosphere to the ocean (or the wind stress) plays an important role in coupled atmosphere-ocean systems. The wind stress t wind is normally parameterized by the neutral 10-m wind velocity U 10N (corrected for stability) and the drag coefficient (C d ) as t wind 5 r air C d U 10N U 10N in model simulations, where r air is air density and U 10N 5 U 10N | |. The drag coefficient C d itself is often parameterized with U 10N , that is, the wind stress magnitude divided by air density, t wind | |=r air 5 u 2 * , is assumed to be a function of U 10N only, where u * is the friction velocity. Since the neutral wind speed profile is logarithmic inside the constant stress layer (but above the direct surface wave impacts), the drag coefficient can be expressed by the roughness length z o as C d 5 k= ln z=z 0 2 , where z 5 10 m and k 5 0.4 is the von Kármán constant, and the normalized roughness length (or the Charnock coefficient) is defined as Z ch 5 z o g=u 2 * (Charnock 1955), where g is gravitational acceleration.
In low to medium wind speeds (approximately U 10N , 20 m s 21 ), previous studies show that C d increases approximately linearly with wind (e.g., Large and Pond 1981;Donelan et al. 2004;Edson et al. 2013). In tropical cyclone (TC) conditions (approximately U 10N . 25 m s 21 ) C d varies widely, roughly from 1.5 3 10 23 to 4.5 3 10 23 (e.g., Powell et al. 2003;Donelan et al. 2004;Sanford et al. 2011;Bryant and Akbar 2016;Hsu et al. 2019). Although some of this large variability of observed C d is likely due to difficulties in measuring wind stress in extreme high-wind environments, it is also expected that the drag coefficient depends on factors other than U 10N , such as surface wave conditions (sea states).
Some previous studies suggest that the drag coefficient depends on sea states in addition to U 10N (e.g., Toba et al. 1990;Smith et al. 1992;Johnson et al. 1998;Taylor and Yelland 2001;Edson et al. 2013), including TC conditions (e.g., Holthuijsen et al. 2012;Chen et al. 2013;Reichl et al. 2014;Chen et al. 2020). Edson et al. (2013) (and studies cited in their study) suggest that the Charnock coefficient (Z ch ) decreases with the wave age (c p =u * , where c p is the wave phase speed at the wave spectral peak) based on observations under low to medium wind conditions. Toba et al. (1990) show that Z ch increases with the wave age if laboratory observations (with very small wave ages) are included. Some studies (Taylor and Yelland 2001;Edson et al. 2013) suggest that Z ch increases with wave steepness (H s /L p , where H s is the significant wave height and L p is the wavelength at the wave spectral peak) in low to medium wind speeds.
Under TC conditions the sea state dependence of C d has been addressed in very few observational studies. Holthuijsen et al. (2012) suggests that swells misaligned with wind impact the C d estimation based on observations under Hurricane Bonnie. Several modeling studies have simulated the sea state-dependent C d under TCs (Chen et al. 2013;Reichl et al. 2014;Chen et al. 2020). The model results typically show that C d in the front-right quadrant of TCs is similar to or higher than that in the front-left quadrant (in the Northern Denotes content that is immediately available upon publication as open access.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-21- 0246.s1.
Hemisphere), which is very different from the spatial pattern observed by Holthuijsen et al. (2012). Swells misaligned with local wind may also cause the misalignment of the wind stress direction from the wind direction (Chen et al. 2013;Reichl et al. 2014;Hsu et al. 2019;Chen et al. 2020).
Although most studies estimate the drag coefficient in TC conditions using atmospheric measurements, some studies (Hsu et al. 2017(Hsu et al. , 2019;;Sanford et al. 2011) estimate the wind stress using the upper-ocean current observations. Specifically, if the horizontal momentum equations for the ocean currents are integrated vertically from the surface down to a certain depth, where the turbulent stress is sufficiently small, the wind stress is approximately equal to the vertical integration of four (time derivative, Coriolis, nonlinear, and pressure gradient) terms, as discussed in more detail in section 4. Sanford et al. (2011) estimate the wind stress from the ocean observations by assuming that the wind stress is approximately equal to the integral of the two linear terms (time derivative and Coriolis terms) in the momentum equation (i.e., the nonlinear and pressure gradient terms are negligible). Their investigation reports a relatively low C d around 1.5 3 10 23 for 10-m wind speed U 10 . 25 m s 21 . Hsu et al. (2017) have significantly improved this approach of wind stress estimation. They first estimate the wind stress from the two linear terms as in Sanford et al. (2011), using ocean current velocity profiles from three electromagnetic autonomous profiling explorer (EM-APEX) floats deployed under Typhoon Megi during the 2010 ITOP (Impact of Typhoons on the Ocean in the Pacific) program (D'Asaro 2014). They restrict their wind stress estimation to the rightfront and left-front quadrants of the storm (i.e., prior to arrival of the storm center), where the nonlinear and pressure gradient terms are relatively small. Next, they estimate the contributions from these two terms (neglected in the initial estimates) by running the Price-Pinkel-Weller 3D (PWP3D) ocean model (Price et al. 1986), and correct the wind stress estimates. Finally, the drag coefficient is estimated using the wind speed field that has been carefully constrained by interpolating two wind field snapshots that are generated using stepped frequency microwave radiometer (SFMR) surveys and dropsonde measurements, also performed during the ITOP program. Their final estimates are in the range of roughly 2 3 10 23 to 3.5 3 10 23 , which are much higher than the earlier estimates by Sanford et al. (2011) and more consistent with other studies.
In their subsequent study (Hsu et al. 2019), they estimate C d in front of five tropical cyclones (Typhoons Megi and Fanapi and Hurricanes Frances, Gustav, and Ike) by using the same approach as in Hsu et al. (2017), using the EM-APEX observations performed during the 2004 CBLAST (Coupled Boundary Layer Air-Sea Transfer) program (Black et al. 2007) and the ITOP program. They report significant stormdependent variations of the drag coefficient, roughly ranging from 1 3 10 23 to 3.5 3 10 23 . They also suggest that C d is strongly dependent on a new parameter "effective wind duration," which is a function of U 10 , storm translation speed, and position relative to the storm center. They find that the faster the storm moves, the larger C d is, and suggest that the dependence they find may be due to variations in the wave age and wave breaking frequencies under different regimes.
In this study, we estimate the wind stress (drag coefficient) by combining the same ocean current observations under the five TCs (Hsu et al. 2019) and a coupled wave-ocean model. We generate the wind fields in a different manner (see section 2a), and estimate the wind stress by directly comparing the observed ocean currents and simulated ocean currents using the coupled model. The three main objectives of this study are 1) to clarify the effect of wind speed differences among different wind products on C d estimates under TCs, 2) to carefully examine and correct for the contributions from the nonlinear and pressure gradient terms in the wind stress estimates using the model simulations, and 3) to perform a comprehensive analysis of the sea state dependence of wind stress and drag coefficient. We also examine in detail various surface wave impacts on the upper ocean currents and wind stress estimates.
Running the coupled wave-ocean model requires a continuous wind field covering the entire model domain over a time period of about 5 days. Therefore, we may not use the same approach of wind field generation by Hsu et al. (2019), relying on wind surveys that are limited in space/time. Instead, a continuous wind field is formed using the parametric wind model (Chen et al. 2020;Reichl et al. 2016a,b) based on the postseason reanalysis best-track database produced by U.S. Navy Joint Typhoon Warning Center (JTWC) for typhoons, and by the National Oceanic and Atmospheric Administration/National Hurricane Center (NOAA/NHC) for hurricanes (HURDAT). The best-track data file contains storm position, propagation speed and direction, center and environmental pressure, radius of outermost closed isobar, maximum wind speed (V max ), radius of maximum wind (R max ), and radii of 18 m s 21 (R 18 ) and 26 m s 21 (R 26 ) winds in the four quadrants of the storm every 6 h, when available.
In this study we assume that the 10-m wind speed U 10 is the neutral 10-m wind speed U 10N , since the stability effect is small in high-wind conditions U 10 . 25 m s 21 .
For Typhoons Megi and Fanapi, a symmetric azimuthal wind field is first generated by the following modified Rankine vortex:
where n is an empirical decay parameter. To introduce a radial wind component and asymmetry to the wind field, the inflow angle (increasing from 08 to 218 with increasing distance from the storm center) is specified and 70% of the storm translation speed is added to the wind field as in Moon et al. (2003). In the western Pacific, there are fewer in situ observations available than in North Atlantic. The best-track parameters from JTWC are generated mostly based on satellite measurements and are not as accurate as those in North Atlantic where in situ observations are routinely available. Therefore, for both typhoons, V max and R max are adjusted from those in the best-track database. Specifically, we first empirically determine the three parameters (V max , R max , and n) when the SFMR survey is available (three surveys for Megi and two surveys for Fanapi) by minimizing the root-mean-square error (RMSE) between the parametric wind speed U 10 and the observed SFMR wind speed U 10 when U 10 . 25 m s 21 .
Figure 1 shows the SFMR survey locations. (The third survey of Megi is not shown since it does not impact our drag coefficient estimates.) The SFMR wind speeds have been corrected using the dropsonde measurements as discussed in Hsu et al. (2017), before they are used for adjusting V max and R max . Once V max and R max are determined, we then calculate the bias correction factor c V (or c R ), which is the ratio of the empirically determined V max (or R max ) relative to the V max (or R max ) in the best-track database. In between the surveys the three parameters (c V, c R , and n) are interpolated. Before the first survey and after the last survey these parameters are set constant. The time series of V max and R max (both before and after the adjustment) as well as n are shown in Fig. 2. The adjusted V max of Megi and Fanapi increase to 86 and 63 m s 21 from 80 and 54 m s 21 , respectively, and R max decreases to 15 and 22 km from 19 and 28 km, respectively. Figure 3 compares the range of parametric model wind speeds (blue area) and the SFMR wind speeds (black dots), both of them plotted against radius (distance from the storm center), during each survey period. The adjusted parametric wind profiles generally agree well with the SFMR winds. (Note that the parametric wind radial profile varies in time and also depends on the direction relative to the storm propagation direction.)
Figure 4 compares all the SFMR wind speed observations and the corresponding parametric model wind speeds (evaluated at the same time and at the same location). Again, they generally agree well, with the RMSE of 5.8 m s 21 for Megi and 4.6 m s 21 for Fanapi. The high wind observations of Fanapi are overestimated by the parametric model. This is likely because the wind speed spatial distribution of the model is strongly constrained and cannot always be matched with observations. The model requires that the maximum wind speed occurs to the right of the storm track (which is statistically true in real storms), but the maximum wind region was observed in front of the storm during the SFMR surveys of Fanapi. Nevertheless, the mean biases between the parametric wind and the SFMR wind for Megi and Fanapi are 1.77 and 1.93 m s 21 , respectively. This bias is well within the expected uncertainty of SFMR observation (Hsu et al. 2017).
Since the parametric wind model does not provide the background wind field away from the storm, the 10-m wind speed field from the Japanese Meteorological Society Reanalysis product (JRA55; Tsujino et al. 2018) is combined with the parametric wind field. Specifically, the parametric wind is replaced by the JRA55 wind when the latter is larger than the former.
For Atlantic hurricanes the best-track data are more reliable because they are produced by combining many observations including satellite measurements, aircraft surveys, and other in situ observations. Therefore, they can be used without corrections. Furthermore, they provide reliable estimates of radii of 26 and 18 m s 21 (R 26 and R 18 ). When the modified Rankine vortex model is fitted to these parameters, it tends to overestimate the wind speed away from the storm center (Fig. 5). We therefore generate the symmetric azimuthal wind field by combing a Rankine vortex and a vortex with exponential decay: V(r) 5
where V 22 5 22 m s 21 , R 22 5 (R 18 1 R 26 )=2, and (R 18 , R 26 ) are the averages of (R 18 , R 26 ) in four quadrants. The decay parameter n is set so that V 5 V 22 at r 5 R 22 . The time series of V max , R max , R 18 and R 26 for Hurricanes Frances, Gustav, and Ike are shown in Fig. 6. As before, the inflow angle and 70% of the translation speed are added to the wind field and the JRA55 U 10 is used for the background wind. The top panels of Fig. 7 shows examples of the two-dimensional wind fields generated by combining the parametric wind model and JRA55 wind [hereafter URI (University of Rhode Island) wind] for five TCs. The timing of each panel has been chosen to be during our wind stress estimation period. They can be compared with the wind fields used by Hsu et al. (2019) [hereafter APL (Applied Physics Laboratory) wind], which has been generated by interpolating the wind maps generated by the SFMR wind surveys. There are significant differences in the spatial patterns of wind speed. Notice, in particular, that the URI wind field is generally broader and more symmetric and the largest wind speed always appears to the right of the storm track, while the APL wind fields are more complex. Since it is not immediately clear which wind products are more accurate, we will use both wind fields in section 5a to investigate the effect of different wind products on the drag coefficient estimates.
During the passage of Frances and Gustav, the National Data Buoy Center (NDBC) buoys (see left panel in Fig. 8) captured the wind field, which can be compared with the URI and APL wind fields. Since the anemometer height varies from 3.6 to 5 m at these NDBC buoys, the 10-m wind speed is calculated by assuming a logarithmic profile (Shearman and Zelenko 1989). The SFMR surveys and best-track message files are both based on 1-min-averaged winds while NDBC buoys report 10-min-averaged winds. Therefore, a gust factor 0.87 is applied to convert the URI and APL winds to the 10-min-averaged winds for comparison (Harper et al. 2010).
Figures 8a and 8b compare the URI wind and observed wind under Frances and Gustav, respectively. Figure 8c compares the APL wind and observed wind under Gustav. (Note that APL wind fields are available only between SFMR surveys.) Generally, both URI wind and APL wind are in good agreement with the NDBC buoy observations except for the low wind speed range. In this study wind stress estimation is performed only at high wind speeds (U 10 . 25 m s 21 ). Unfortunately, there are not enough measurements during such high wind speeds from the NDBC buoy observations to validate the wind models.
In this study observations from 12 EM-APEX floats deployed under five TCs are used, including three floats under Typhoon Megi, four floats under Typhoon Fanapi, three floats under Hurricane Frances, and two floats under Hurricane Gustav (Table 1 and Fig. 1). One float (float 3766) deployed on the left side of Gustav's track drifted to the right side of Ike's track and provided the measurements under Ike as well (Figs. 1 and 7). These EM-APEX float observations provide the vertical profiles of horizontal current, temperature, and salinity roughly every 30 min (Sanford et al. 2011;Hsu et al. FIG. 4. Comparison of URI U 10 and SFMR observed U 10 for all SFMR surveys under Typhoons (left) Megi and (right) Fanapi. Z H O U E T A L .
2017, 2019). The deployment time of each EM-APEX float was about 1 day before the storm arrival (Table 1). The detailed trajectory of 12 EM-APEX floats under TCs are shown in Fig. 7.
As discussed by Hsu et al. (2017Hsu et al. ( , 2019)), the measured ocean currents include wind driven currents as well as tidal currents and low-frequency background currents. Since we estimate the wind stress by comparing the observed currents to model simulated currents that are purely wind driven, both the tidal currents and background currents are removed from the measured currents. The full description of extracting wind driven current velocity from the direct EM-APEX velocity measurements is given in Hsu et al. (2017Hsu et al. ( , 2019)). Following Hsu et al. (2017Hsu et al. ( , 2019)), missing measurements in the upper 30 m under Megi and Frances are filled with the uppermost velocity measurements.
The coupled Modular Ocean Model 6-WAVEWATCH III (MOM6-WW3) model is used in this study to simulate the coupled wave-ocean system under TC winds. The MOM6 is FIG. 6. Time series of (left) maximum wind speed, (middle left) radius of maximum wind, and radii of (middle right) 26 and (right) 18 m s 21 wind speeds for Hurricanes (top) Frances, (middle) Gustav, and (bottom) Ike every 6 h. Two vertical gray lines indicate the time period of C d estimation in this study.
the latest version of the modular ocean model with a hybrid vertical coordinate system developed and maintained by NOAA Geophysical Fluid Dynamics Laboratory (Adcroft et al. 2019). For the wave simulation the WW3 version 6.07 is used (WW3DG 2019). The wind input term (S in ) and white-capping terms (S dis ) in WW3 in this study are setup following the WW3 ST4 version in Liu et al. (2017) and Chen et al. (2020), which showed good agreement between model results and observations under Hurricane Ivan (see Supporting Information D in the online supplemental material for more details). In the coupled MOM6-WW3 system, the WW3 provides the dominate wavelength and the Stokes drift vector (u S ) to the ocean model. The MOM6 provides the near surface horizontal current vector (u W ) to the WW3. In these simulations we use the K-profile parameterization (KPP) for ocean surface boundary layer mixing FIG. 8. (left) Locations of NDBC buoys (black triangles) in North Atlantic and Gulf of Mexico and tracks of Frances (yellow dotted curves) and Gustav (purple dotted curves). Comparison of URI U 10 and buoy observations under (a) Frances and (b) Gustav. (c) Comparison of APL U 10 and buoy observations under Gustav. Blue symbols in (b) are DATA when APL wind was not available. Comparison of significant wave height (H s ) between simulation and observation under (d) Frances and (e) Gustav. In all panels different symbols indicate different buoys. The root-mean-square error and bias between model simulation and observation are labeled.
in MOM6 via the public Community Vertical Mixing (CVMix) project.
Three large, deep-water (4000 m) computational domains are used for both the ocean and the wave models. As shown in Fig. 1, the domain in the western Pacific for Megi and Fanapi is 108-408N, 1108-1408E. The domain in the North Atlantic for Frances is 158-308N, 908-558W. In Gulf of Mexico, the domain for both Gustav and Ike is 158-308N, 1008-658W. The spatial resolution is 1/248, which is around 4.5 km, and the temporal resolution is 300 s for both MOM6 and WW3. We note that it is important to use a spatial resolution of 1/248 or higher for WW3, since Chen et al. (2018) found that a coarser spatial resolution of WW3 may introduce a significant bias of wave spectra under TCs. The vertical resolution in MOM6 is 2 m in the upper ocean (50 m), increases linearly to 10 m at a 100-m depth and gradually increases below. The vertical coordinate in MOM6 is set to the stretched geopotential, or Z * mode for these experiments. The surface wave spectrum in WW3 is discretized using 48 directions and 40 relative frequencies. The initial temperature and salinity profiles in the model experiments are prescribed spatially homogeneous, and are determined by spatially and temporally averaging the vertical profiles from EM-APEX floats observations before the arrival of TC's eye, when JRA55 U 10 is larger than the parametric U 10 .
Since our later investigation of the sea state dependence of drag coefficient (section 6) relies on the WW3 simulation results, it is important to clarify their accuracy. The wave parameters we will later utilize include misalignment angle c between wind and dominant (or peak) waves, wave age c p =u * , and wave steepness H s /L p , in conditions of dominant wavelength from 129 to 384 m and H s from 6.5 to 15 m. Therefore, WW3 predictions of the dominant wave direction, dominant wavelength, and significant wave height need to be sufficiently accurate in such conditions. In three previous studies (Moon et al. 2003;Fan et al. 2009;Liu et al. 2017) WW3 predictions of these wave parameters under TCs were carefully validated against National Aeronautics and Space Administration (NASA) Scanning Radar Altimeter (SRA) measurements. The accuracy of wave observations of SRA has been well established (Walsh et al. 1985(Walsh et al. , 1987(Walsh et al. , 1989). In the 1990s, the SRA mode of the NASA 36-GHz multimode airborne radar altimeter was developed for installation on NOAA research aircraft that were flying missions into tropical cyclones (PopStefanija et al. 2021). The detailed data processing of SRA measurements are well documented in Walsh et al. (2002) and Wright et al. (2001). In particular, SRA observations of dominant wavelength and direction are highly reliable provided the dominant waves are longer than 100 m (E. Walsh 2022, personal communication). Moon et al. (2003) conducted detailed comparison of wave spectra from WW3 with SRA observation under Hurricane Bonnie (1998). Their Fig. 13 (reproduced in the supplemental material, Supporting Information A) shows excellent agreement of dominant wavelength, dominant wave direction and significant wave height along all SRA tracks in open ocean, except for one SRA track near the landfalling storm. Using SRA measurements under Hurricane Ivan (2004), Fan et al. (2009) also show excellent agreement of dominant wavelength, dominant wave direction and significant wave height along SRA tracks in their Fig. 4 (reproduced in the supplemental material, Supporting Information A). In particular, they show that the agreement of significant wave height improves when ocean surface currents are included in the WW3 simulation. The coupled wave-ocean model in this study also includes the ocean current effect on wave simulations. Fan et al. (2009) also show good agreement of significant wave height (above ∼6 m) between SRA measurements and satellite measurements from Envisat and ERS-2. More recently, Liu et al. (2017) investigated WW3 performance under Hurricane Ivan ( 2004) by comparing WW3 with updated source terms with SRA observations. Their Fig. 5 clearly confirms excellent agreement of significant wave height, mean wave direction, and mean wave period.
Based on the comparison between WW3 simulations and SRA measurements in these three studies (Moon et al. 2003;Fan et al. 2009;Liu et al. 2017), we may estimate the WW3 uncertainty of dominant wave direction to be ∼108. The uncertainty of dominant wavelength is ∼50 m, which introduces ∼0.1 uncertainty of wave age. For significant wave height above 6 m, its uncertainty of ∼1 m introduces ∼0.005 uncertainty of wave steepness. These uncertainty estimates are sufficiently small for our later analysis of the sea state dependence of drag coefficient in section 6.
Finally, Figs. 8e and 8d confirm that the WW3 simulations of significant wave height under high-wind condition are consistent with the buoy observations under Hurricane Gustav. The bias between WW3 simulated H s and NDBC buoy observations are comparable to those in previous studies (Fan et al. 2009;Liu et al. 2017;Moon et al. 2003).
In this study we investigate the impacts of ocean surface waves (sea states) on four different air-sea interaction processes. First, surface waves modify the drag coefficient (or the equivalent surface roughness) mainly because wind stress is mostly determined by form drag of a spectrum of surface waves in higher wind speeds. Second, when surface waves are growing/decaying, the stress applied to the upper ocean may be different from the wind stress (wave-induced momentum flux budget). Third, upper-ocean turbulent mixing is modified by surface waves (Langmuir turbulence). Fourth, upper-ocean currents interact with surface waves through Coriolis-Stokes force, Stokes-shear force, and other processes (Suzuki and Fox-Kemper 2016). While the first wave effect (sea statedependent drag coefficient) affects both atmosphere and ocean, the other three wave effects only impact the subsurface currents and turbulence.
In this study we thoroughly investigate the first wave effect (sea state dependence of drag coefficient) in section 6. While the other three wave effects do not impact the drag coefficient itself, they may impact our drag coefficient "estimation," because our estimation is made using the upper-ocean current observations. Therefore, a series of numerical experiments are performed to examine the other three surface wave impacts, that is, Langmuir turbulence, the wave-induced momentum flux budget, and wave-current interactions, on the upper-ocean current responses and the wind stress estimation (Table 2). In MOM6 the following Boussinesq momentum equations (given here in geopotential/height coordinates) are solved (Adcroft et al. 2019):
where u is the (3D) Eulerian velocity, f is the local Coriolis frequency, r 0 is the constant Boussinesq reference density, and r is the in situ density. In the horizontal momentum equations F represents the accelerations due to the divergence of stresses. In this study the horizontal stress is negligible and F is practically equal to the vertical gradient of the vertical frictional turbulent stress t i /z. The turbulent stress t i is parameterized using two different KPP schemes (KPP-iLT and KPP-LT) proposed by Reichl et al. (2016b), who have carefully tuned the KPP schemes against a large number of large-eddy simulation (LES) runs under TCs. The KPP-iLT is a standard KPP model without explicit surface wave impacts, but the critical Richardson number is set Ri cr 5 0.35 so that it is consistent with typical (sea state independent) Langmuir turbulence conditions under TCs. The KPP-LT incorporates explicit sea state-dependent Langmuir turbulence, mainly by setting the eddy viscosity as a function of the turbulent Langmuir number (ratio of the Stokes drift velocity to the friction velocity), and by parameterizing the turbulent mixing in terms of the Lagrangian current shear instead of the Eulerian current shear. Further details of these mixing schemes are given in Reichl et al. (2016b).
As shown in Table 2, Experiment A is performed with the implicit (sea state independent) Langmuir turbulence KPP model (KPP-iLT), while Experiment B is performed with the explicit (sea state dependent) Langmuir turbulence KPP model (KPP-LT). Experiment C is carried out to study the influence of the air-sea momentum flux budget, that is, the difference between the momentum flux from atmosphere (wind stress) and the momentum flux into upper ocean, due to growing/decaying surface waves. In this experiment the surface forcing for the upper ocean is modified from the wind stress by Dt bg i (Fan et al. 2009(Fan et al. , 2010;;Smith 2006):
Here, M i 5 0 2' u S i (z)dz is the vertically integrated Stokes drift, S ij is the radiation stress, and u w i is the near surface current (Smith 2006). Fan et al. (2010) investigated the air-sea momentum flux budget by calculating the first two terms on the right of Eq. ( 5) since the last two terms are much smallerroughly by the ratio of the surface current and the group velocity of dominant waves. Nevertheless, these smaller terms are retained in this study because they are comparable to the horizontal Stokes shear force discussed below. The mixing scheme in Experiment C is KPP-LT.
Experiment D includes the wave-current interaction effects in addition to the KPP-LT and the air-sea momentum flux budget, that is, includes all the wave impacts addressed in this study. Namely, the governing equation for the ocean currents are modified to
which include three surface wave effects (Suzuki and Fox-Kemper 2016). Namely, 1) the advection is by the Lagrangian current (u L 5 u 1 u S ) instead of by the Eulerian current [second term in Eq. ( 6)], 2) the Coriolis force applies to the Lagrangian current instead of the Eulerian current [Coriolis-Stokes force, third term in Eq. ( 6)], and 3) the Stokes shear force applies [right-hand side terms in Eq. ( 6) and ( 7)]. The continuity equation is also rewritten in terms of Lagrangian current instead of Eulerian current.
In principle, the MOM6 code needs to be modified to properly solve these equations. However, in this study we take a simpler (but less accurate) approach. First, Eq. ( 6) is rewritten in terms of the Lagrangian current:
Then, the left terms of Eqs. ( 8) and ( 7) are identical to the momentum equations solved in MOM6 [Eqs.
(3) and ( 4)], if the MOM6 solves for the Lagrangian current u L instead of the Eulerian current u, and F is also parameterized using u L . For the horizontal momentum equation [Eq. ( 8)] only the three extra terms on the right need to be added as extra forcing. Since the Stokes drift decays with depth rather quickly, these terms only apply to the shallow upper layer. We therefore assume that these three forcing terms apply at the surface instead of in the interior. Specifically, if we integrate Eq. ( 8) from depth 2H s to the sea surface, where H s is the depth of the Stokes layer (where the Stokes drift is significant),
The first term (integral of three terms) on the right of Eq. ( 9) is denoted as Dt wc i and treated as a modification to the wind stress. The vertical momentum equation, Eq. ( 4), is also modified by the vertical Stokes shear force as shown in Eq. ( 7). The vertically integrated Eq. ( 7) becomes
that is, the pressure field below the Stokes layer is modified as if the surface atmospheric pressure is modified by the last term on the right of Eq. ( 10), which is named Dp wc . This surface pressure correction term Dp wc is also introduce into Experiment D. Although our approach does not accurately resolve the surface wave impact inside the Stokes layer, we consider this as a reasonable approach (as a first step) to perform a sensitivity study of the wave-current interaction effects.
In this section we describe how we estimate the wind stress by combining the observation and the model.
As discussed by Sanford et al. (2011) and Hsu et al. (2017Hsu et al. ( , 2019)), the depth-integrated horizontal momentum equation from the sea surface to a depth 2H becomes
Here, t iz50 and t iz52H are the surface wind stress (or the modified stress if the surface wave effects are included) and the stress at z 5 2H. If we choose a sufficiently deep 2H so that the stress is negligible there, the wind stress is equal to the vertical integral of four terms on the left hand side of Eq. ( 11). In this study we choose H 5 100 m as in Hsu et al. (2019). Hsu et al. (2017Hsu et al. ( , 2019) ) first estimate the wind stress by integrating the two linear terms ( t u i 1 « i3j fu j ) using the EM-APEX ocean current observations. Since the nonlinear and pressure gradient terms cannot be obtained from the observational data, they use model-simulated nonlinear and pressure gradient terms to estimate their contribution to the stress calculation, and restrict their wind stress estimation to the front of TCs, where these terms are small. In Hsu et al. (2017Hsu et al. ( , 2019) ) the wind stress is estimated for each EM-APEX profile, but the time derivative t u i is estimated from two profiles (one before and one after) that are Dt apart. Since the float drifts over that time period, the measurements actually provide Du i /Dt, where Du i is the difference of the observed horizontal current (including the float advection effect), and it is different from t u i . We investigate the difference between Du i /Dt and t u i in section 5.
Following Hsu et al. (2017Hsu et al. ( , 2019)), we first estimate the wind stress r obs 5 (s obs
x , s obs y ) directly from the wind-driven current observations without correcting for the nonlinear and pressure gradient terms: Here, (u 10 , y 10 ) are the zonal and meridional components of U 10 , and an asterisk denotes the uncorrected estimates. In section 5a we perform this exercise using both the URI and APL wind fields in front of TCs to investigate the effect of wind speed uncertainty on the (uncorrected) drag coefficient estimates.
In this study the effects of nonlinear and pressure gradient terms are not introduced as corrections to be added. Instead, we will take the following steps.
1) The coupled MOM6-WW3 model is run using the URI wind field and an assumed (initial guess) drag coefficient (C init d ) as a function of U 10 . 2) We simulate the (uncorrected) stress calculation from the float observations by using the model simulated current at the same time and at the same location as the float observation. Specifically, we calculate r sim 5 (s sim x , s sim y ) with
for each EM-APEX profile, where u sim 5 u sim , y sim ( ) is the simulated current.
3) We compare the magnitude and direction of r obs and r sim for each data point (i.e., each EM-APEX profile), including those behind the storm. Specifically, we calculate the magnitude ratio M r :
and the misalignment angle f between the simulation and the observation,
f 5 tan 21 s obs y s obs x 2 tan 21 s sim y s sim x :
Assuming that the URI wind fields are reasonably accurate, if the drag coefficient C init d used in the simulation is correct, we would expect that M r is close to 1 and f is close to 0 on average. (If not, we need to modify C init d and repeat the simulations. However, we will find that our initial assumed drag coefficient is sufficiently accurate.) Furthermore, if we assume that r sim responds linearly to the locally imposed wind stress (i.e., |r sim | is proportional to the magnitude of the imposed stress, and the direction of r sim changes in the same manner as the direction of the imposed stress if it is changed), the drag coefficient can be estimated from each observation as
f is the estimated misalignment angle between the wind speed and the wind stress, and the downwind drag coefficient is calculated as C d|| 5 C d cosf (Hsu et al. 2017). (If f is positive, the wind stress direction is to the left of the wind speed direction.) Then, the friction velocity u * , the surface roughness z o , and the Charnock coefficient Z ch are calculated using u 2 * 5 C d U 2 10 , C d 5 k= ln z=z 0 2 with z 5 10 m, and Z ch 5 z o g=u 2 * . As an initial estimate of the drag coefficient, we use the drag coefficient implemented in the Geophysical Fluid Dynamics Laboratory (GFDL) hurricane model in 2015 (Ginis et al. 2015) (hereafter GFDL C d ; Fig. 9). This particular C d parameterization is chosen because it produced the best tropical cyclone intensity forecasts using the GFDL and HWRF hurricane models based on a large number of storm simulations (Ginis et al. 2015;Biswas 2018). It is close to the COARE 3.5 formulation at wind speeds lower than 23 m s 21 (Edson et al. 2013) and consistent with the theoretical formulation of Soloviev et al. (2014) above 30 m s 21 . Figure 9 shows that the GFDL C d is close to the average of the previous observations under TCs including the results of Hsu et al. (2019). We find that the model simulation using this drag coefficient is reasonably accurate, i.e., the resulting bin-averaged values of M r and f are sufficiently close to 1 and 0, respectively. Therefore, we do not modify C init simulation. However, in one experiment we have estimated C d based on ocean simulations using a different C init d (1.2 times GFDL C d ). We have found that the resulting C d estimates are almost identical to the original estimates, that is, our C d estimation is not sensitive to the choice of C init d (see supplemental material, Supporting Information B).
Our approach of the drag coefficient estimation described in the previous subsection relies on an assumption that r sim responds linearly to the locally imposed wind stress, that is, r sim is close to the imposed wind stress and the contribution from the other terms is reasonably small. Therefore, we carefully examine the validity of this assumption for each observational data point (each EM-APEX profile). Namely, 1) We repeat all the simulations with a drag coefficient that is 20% larger. If the resulting magnitude r sim increases more than 30% or less than 10% for a particular data point, it is excluded from the analysis. 2) If r sim is less than 70% or more than 130% of the imposed wind stress magnitude for a particular data point, it is excluded from the analysis.
Note that these two steps are purely based on the model simulations; we do not impose any conditions on the observational data. After these quality control processes, the number of data samples with U 10 $ 25 m s 21 has decreased from 326 to 195. Nevertheless, we now utilize a significant number of data points behind the storm that were not included in the previous studies.
By repeating the analyses using increased and reduced data samples (with more relaxed and more strict data quality control criteria), we have ascertained that adjusting these criteria do not significantly change the conclusions of this study (see supplemental material, Supporting Information C).
In this section, the results are mostly presented as bin averages for wind bins of 25-30, 30-35, 35-40, 40-45, and .45 m s 21 (from Fig. 10). The horizontal error bar is the standard deviation of wind speed in each bin. The vertical error bar is the 95% confidence level in each bin. The 95% confidence level is infinite when there is only one data point within the bin.
The accuracy of estimated drag coefficient depends on the accuracy of wind speed products. To investigate the impact of different wind products on drag coefficient estimates, we present C * d , C * d|| 5 C * d cosf * , and f * with both the URI wind and APL wind in Fig. 10. Note that C * d , C * d|| , and f * are estimated from the integration of linear (measurable) terms only and are not as accurate as the estimates presented later. The FIG. 11. As in Fig. 10, but separated for each TC and in front-right sector. There are 109 data points in the front-right sector and 30 data points in the front-left sector under all five TCs, since we have not applied our data quality criteria for this exercise.
In general, the results of C * d , C * d|| , and f * from the two wind products agree reasonably well in both the front-right and front-left sectors (Fig. 10). The vertical error bars (95% confidence) significantly overlap for most data points. This suggests that our drag coefficient estimates are not significantly affected by wind speed uncertainty if a large number of data points are averaged. Notice that both C * d and C * d|| are significantly lower in the front-left sector than in the front-right sector. This is likely because dominant wave direction is more misaligned from wind direction in the front-left sector, as discussed in section 3c. (In front-left and front-right of TCs, roughly 80% and 30% of data show large misalignment angles exceeding 458, respectively.)
If the results are shown separately for each TC, the difference between the two wind products becomes more noticeable (Figs. 11 and 12). In particular, significant differences appear in lower wind conditions (|U 10 | , 30 m s 21 ) under Megi and Frances, corresponding to the noticeable differences between the URI and APL wind fields (Fig. 7). The figures also confirm that C * d and C * d|| significantly vary from storm to storm, as pointed out by Hsu et al. (2019). We suspect that the wind speed differences among different wind products is partly responsible for the observed storm dependent variations of C * d and C * d|| .
An example of observed and simulated wind driven current time series is shown in Fig. 13. We have chosen Float 3763 under Typhoon Megi since it passed under high wind speeds (Fig. 14a). The model simulations generally agree well with the observations. The initial excitation of wind-driven current and its rotation in later times are well reproduced in the simulations. The differences among the four model results (Experiments A-D, discussed in section 3) are subtle, i.e., the surface wave impacts are not large. Note that both the KPP-iLT (Experiment A) and KPP-LT (Experiments B-D) have been tuned against the same set of LES runs including the Langmuir turbulence effects (Reichl et al. 2016b). If the simulation is performed without the enhanced mixing by Langmuir turbulence, the difference is expected to be more significant (Reichl et al. 2016a). The currents in Experiment D (with the wave-current interactions included) appear to be slightly weaker. We will later find that our drag coefficient estimates are slightly larger with Experiment D as a result.
Next, we investigate the results of the data quality control process discussed in section 4c. Five examples (one example from each TC) of the time series of vertically integrated momentum budget terms are presented in Fig. 14. In this figure the dashed cyan, dashed purple, and dashed green lines show the magnitudes of combined linear (first and third) terms, nonlinear (second) term, and pressure gradient (fourth) term, on the right side of Eq. ( 11), respectively. If these four terms (as vectors) are added, it is equal to the wind stress vector [Eq. ( 11)], whose magnitude is shown by the solid black line. In addition, the time series of r obs and r sim are shown by the solid gray and solid red lines, respectively. Recall that one of our imposed criteria is that r sim (solid red line) is between 70% and 130% of the imposed wind stress magnitude (solid black line). The figure clearly shows that our data quality control process has successfully removed data points when the nonlinear (dashed purple) and/or pressure gradient (dashed green) terms are significant. Even with this strict data control a significant number (195) of data points have been found usable, including many in the rear of storms (positive arrival time in Fig. 14), which are excluded in Hsu et al. (2017Hsu et al. ( , 2019)). It is noteworthy that significant differences sometimes appear between the integrated linear terms (dashed cyan) and r sim , that is, the time derivative calculated from the two float profiles can be significantly different from the partial time derivative because of the drifting of the float.
We now present the results of estimated M r , C d , C d|| , and f using all 195 data points under five TCs in Fig. 15 (top panels). The magnitude ratio M r between the observed r obs and the simulated r sim (first panel from left) is reasonably close to 1 except for the highest wind bin of U 10 . 45 m s 21 . The number of data points in this bin is small and the error bar is large. (The lower end of the error bar is less than 1.) Therefore, it is not clear whether the difference between r obs and r sim is statistically significant. Since the result of M r suggests that the GFDL C d is reasonably consistent with the observation, we do not repeat the model simulation with a modified C init d . Instead, the drag coefficient is estimated by multiplying the GFDL C d and the observed M r , and is shown in the second panel from the left in Fig. 15. The results confirm that our estimated C d is close to FIG. 14. Time series of momentum budget terms. Black solid line shows imposed wind stress magnitude. Cyan, purple, and green dashed lines show combined linear, nonlinear, and pressure gradient terms (all vertically integrated), respectively. Gray and red solid lines show r obs and r sim , respectively. Circle symbols indicate data used in this study. Cross symbols indicate data removed by the data quality control process. Float name is labeled in each panel. FIG. 15. Bin-averaged estimates of (a) M r , (b) drag coefficient C d , (c) along-wind drag coefficient C d|| , and (d) misalignment angle f between wind stress and wind speed, from five TCs. Different symbols represent results from four model experiments. Black lines indicate GFDL C d . Horizontal error bars represent the standard deviations of U 10 . Vertical error bars indicate 95% confidence level. (e)-(h) As in (a)-(d), but exclude data with significantly misaligned swell (|c| . 458).
the GFDL C d . Both our C d and the GFDL C d show weak wind speed dependence and their values are mostly in the range of 2-3 3 10 23 for wind speeds 25-55 m s 21 . The reasonable agreement between our C d estimates and the GFDL C d is particularly encouraging because these two have been derived using completely different approaches; the former is based on in situ upper-ocean observations and the latter is based on a larger number of atmospheric TC simulations (Ginis et al. 2015;Biswas 2018). The observed weak wind speed dependence of our C d (slight decrease first and then increase with U 10 ) appears to be opposite of that of the GFDL C d (slight increase first and then decrease with U 10 ), and this C d trend appears significant because the error bars due to scatter of observational data are relatively small. However, the observed weak C d trend is likely insignificant because using different wind products can affect the C d results by 0.5 3 10 23 or more as discussed earlier (also see Fig. 6 in the supplemental material). In addition, we have repeated C d estimations using more strict and less strict data quality control criteria. While the overall C d values remain very similar (2-3 3 10 23 ), the C d trend can vary depending on the criteria (see supplemental material, Supporting Information C).
The misalignment angle f between r obs and r sim is shown in the rightmost panel in Fig. 15. This is an estimate of the misalignment of wind stress direction relative to wind speed direction. Generally, the estimated value is not large. The error bars mostly intersect f 5 08. However, f tends to be negative (wind stress direction is to the right of wind speed direction) in lower wind speeds. We will later find that the negative f may be related to the negative c (dominant wave direction is to the right of wind speed direction), where c is the misalignment angle of dominant surface wave direction relative to wind speed direction.
The difference of C d estimates using the four different experiments (A-D) is mostly negligibly small. Including different surface wave effects does not appear to influence the C d estimates. The only exception is the highest wind speed bin, where the result of Experiment D is significantly larger than the rest, which is consistent with the earlier observation that the wind driven currents are weaker in Experiment D. This suggests that the wave-current interaction effect can be important in some conditions and that it is desirable to investigate the effect more accurately by modifying the ocean model momentum equations directly, instead of approximating the effect in the boundary conditions as done in this study.
As discussed earlier, our C d estimates (around 2.7 3 10 23 ) in the highest wind bin (U 10 . 45 m s 21 ) is larger than the GFDL C d (around 2.1 3 10 23 ). Furthermore, the estimates by Hsu et al. (2017), obtained using the same observational data, are around 2 3 10 23 and are also lower than our estimates. Therefore, we have carefully investigated the possible causes of this difference between our estimate and the estimates by Hsu et al. (2017). The first cause is the difference in the wind fields. Figure 10a shows that the estimates with the URI wind is larger than those with the APL wind by about 0.3 3 10 23 . Another cause is the different approaches taken to correct for the nonlinear and pressure gradient terms. Hsu et al. (2017) estimates that this correction reduces the drag coefficient estimates, but our momentum budget analysis under Megi (Fig. 14a) shows that the correction increases the estimates (from solid red line to solid black line) in high-wind conditions. Interestingly, our analyses under Fanapi (Fig. 14b) and Frances (Fig. 14c) show that the correction can be opposite, i.e., can decrease the C d estimates as in Hsu et al. (2017). (In our analysis of the highest wind bin, four data points come from Megi, two from Frances, and one from Fanapi.) These analyses suggest that both the wind field specification and the C d estimate correction are difficult to constrain, and that the difference between our C d estimates and the GFDL C d or the C d estimates by Hsu et al. (2017) are not conclusive.
Next, the results of estimated M r , C d , C d|| , and f are presented separately for each TC in Fig. 16. Although the error bars mostly intersect the GFDL C d , the estimates appear to vary from storm to storm. In particular, the bin averages under Gustav are significantly below the GFDL C d . As discussed in section 3a we suspect that this storm dependence is partly due to the wind field uncertainty. It is also possible that some variations are caused by different sea states as discussed in the next section.
Some bin averages show more noticeable differences among the four experiments (A-D), compared to the results combining all TCs in Fig. 15. This suggests that the surface wave effects may impact the upper-ocean responses and the C d estimates in some cases, but their impacts tend to cancel out if a large number of observations are averaged.
In this section we only show the analysis based on the simulation from Experiment C, since the sea state dependence analysis for all 4 experiment are almost identical.
As discussed in section 1 many previous studies suggest that C d is not a function of U 10 only, but depends on sea states. In particular, in TC conditions it is common to observe dominant surface waves that propagate in a different direction from the local wind direction (misaligned swell). Both observations (Holthuijsen et al. 2012) and modeling studies (Chen et al. 2013;Reichl et al. 2014;Chen et al. 2020) suggest that the drag coefficient may be strongly influenced by misaligned swell. Note that the term "swell" is used here for waves that are not locally generated. In low to medium wind speeds the wind forcing on swells is weak or even negative. However, swells can be strongly forced by wind in TC conditions if they are aligned with wind.
We first examine the misalignment angle c between the dominant wave direction and the wind direction for all the data points. In Fig. 17a the spatial distribution of |c| relative to the storm center is shown in the normalized coordinate system, in which the distance from the storm center is normalized by the radius of maximum wind of the TC. As expected, the misalignment is mostly small to the right and right rear of the storm, but is large in the far front (approximately .2-3R max ) and in the left-front quadrant, often exceeding 458.
In Fig. 18 the dependence of C d , Z ch , and f on the misalignment angle |c| is presented. The data are bin averaged every 22.58 increment of |c| for each wind speed bin (color coded filled circles). For C d and f the bin averages are also shown for all wind speeds combined (black empty squares). The two data points with |c| . 908 are excluded. One immediately notices that both C d and Z ch rapidly decrease as |c| exceeds around 458. This strongly suggests that dominant waves misaligned by more than 458 from wind have a negative impact on the drag coefficient. The existing models of sea statedependent drag coefficient (Chen et al. 2013;Reichl et al. 2014;Chen et al. 2020) also predict slight reduction of C d as the misalignment angle |c| increases because the form drag of swell decreases. However, the observed reduction of C d is significantly stronger than the model predictions.
With smaller misalignment between dominant waves and wind speed (|c| , 458), the misalignment f between the wind stress and wind speed scatter above and below 08, but f becomes more consistently negative for |c| . 458. We find that significantly misaligned dominant waves (|c| . 458) almost always propagate to the right of wind (c , 2458). Therefore, there appears to be some correlation between the negative c and negative f, that is, misaligned swell propagating to the right of wind tend to turn the wind stress direction to the right toward the swell direction. This is qualitatively consistent with previous modeling studies (Chen et al. 2013;Reichl et al. 2014;Chen et al. 2020).
As discussed earlier, large misalignment |c| . 458 mainly occurs in the far front and in the left-front quadrant (Fig. 17a), where the drag coefficient appears to be reduced according to Fig. 18a. This particular spatial dependence of |c| is a possible reason why C d in the previous study (Hsu et al. 2019) and our C * d estimates (section 3a) show significant spatial dependence, i.e., they are lower in the far front and in the left-front quadrant.
Since large misalignment |c| . 458 is common under TCs, it is likely that misaligned swells contribute to overall reduction of the mean drag coefficient averaged over all sea states. In to the results including misaligned swell (top panels). This suggests that common occurrence of misaligned swell may be partially responsible for the overall reduction of the drag coefficient under TCs. We also find that the observations under FIG. 18. Bin-averaged estimates of (a) C d , (b) Charnock coefficient Z ch , and (c) misalignment angle f between wind stress and wind speed, plotted against misalignment angle |c| between dominate waves and wind speed.
Hurricane Gustav contain many data points with large misalignment |c| . 458 in lower winds, and they likely contribute to the reduced C d (Fig. 16).
The common parameters to distinguish different sea states in previous studies include wave steepness H s /L p , wave age c p =u * , and input wave age c pi =u * , where c pi is the input phase speed at the input peak frequency that is determined by wind energy input calculated in WW3. Unlike c p , c pi excludes the contribution from waves that are not wind forced (e.g., swell in low to medium winds, significantly misaligned swell in TCs). Therefore, c pi is always smaller than c p but it becomes close to c p in pure wind seas.
In this subsection we investigate the dependence of the estimated C d , Z ch , and f on these sea state parameters. To avoid the self-correlation due to the large scatter of u * , we redefine the wave age and input wave age as c p /U 10 , and c pi /U 10 , respectively.
The spatial distribution of c p /U 10 , c pi /U 10 , and H s /L p are shown in Figs. 17b-d, respectively. In the plots of c p /U 10 and H s /L p distributions, data points with large misalignment |c| . 458 are marked by crosses instead of circles, and are excluded in the analyses of c p /U 10 dependence and H s /L p dependence. However, we retain these data in the analysis of c pi /U 10 dependence because the input wave age should be able to exclude the contamination from misaligned swell.
Figure 17b shows that larger wave ages are more common from the front left to the right of the storm, frequently with misaligned dominant waves, while young seas appear just behind the storm. The distribution of the input wave age in Fig. 17c shows that the wind sea part of the wave spectrum (locally wind forced waves) is mostly young (input wave age less than 0.6). The wave steepness H s /L p appears to be inversely correlated with the wave age c p /U 10 (cf. Figs. 17b,d), and the steepest waves appear just behind the storm. In Fig. 19 the observed dependence of C d , Z ch , and f on the wave age c p /U 10 is presented. Recall that the data points with significantly misaligned swell are excluded. The data are bin averaged every 0.15 increment of c p /U 10 for each wind speed bin (color coded filled circles). For C d and f the bin averages are also shown for all wind speeds combined (black empty squares). As expected, the wave age itself is strongly correlated with wind speed (Fig. 20a). Seas tend to be younger in higher wind speeds. Therefore, it is not easy to distinguish the dependence on c p /U 10 from the dependence on U 10 . Nevertheless, for older seas (c p /U 10 . 0.6) both C d and Z ch appear to decrease with the wave age c p /U 10 , even within the same wind bins. In particular, the dependence of Z ch on c p /U 10 at the lowest wind bin (blue) appear quite similar to the dependence observed in previous studies in low to medium wind conditions (Smith et al. 1992;Johnson et al. 1998;Oost et al. 2002;Edson et al. 2013). Note that the parameterizations proposed in the previous studies based on c p =u * have been converted to those based on c p /U 10 in Fig. 19b. The misalignment angle f becomes significantly negative with increasing wave age for older seas (c p /U 10 . 0.5) (Fig. 19c). Although we have removed the data points with significantly misaligned swell (|c| . 458) in this analysis, older seas often contain dominant waves propagating to the right of wind (misaligned by less than 458) and may contribute to the turning of the wind stress direction to the right.
Interestingly, the dependence of C d and Z ch on the wave age c p /U 10 seems to disappear or even be reversed, i.e., they may increase with increasing c p /U 10 if seas are younger (c p /U 10 , 0.6). Such young waves are rarely observed in low to medium wind speeds; they are observed only with very short fetch, such as in laboratory wind wave flumes. However, they are quite common in TC winds. Although this increasing dependence of Z ch with increasing c p /U 10 is qualitatively consistent with the parameterization by Toba et al. (1990), which is based on combined field and laboratory observations in low to medium wind speeds, our observed Z ch values are significantly lower. This suggests that the high-wind TC conditions and strongly forced laboratory conditions are very different regimes (the same parameterization does not apply).
We next investigate the observed dependence of C d , Z ch , and f on the input wave age c pi /U 10 in Fig. 21. Recall that we have now included the data points with significantly misaligned swell. Since c pi /U 10 is always smaller than c p /U 10 , the overall range of c pi /U 10 is significantly narrower than that of c p /U 10 . Nevertheless, we still observe strong correlation between c pi /U 10 and U 10 (Fig. 20b); the wind sea part of wave spectrum is younger with higher wind speed. As in the case of wave age dependence, C d and Z ch seem to remain constant or slightly increase with the input wave age c pi /U 10 for younger seas. However, their overall dependence on c pi /U 10 is weak, that is, c pi /U 10 does not show clear sea state dependence in this wind speed range.
Last, the dependence of C d , Z ch , and f on the wave steepness H s /L p is shown in Fig. 22. Here, the data with significantly misaligned swell have been excluded. The steepness H s /L p is also correlated with wind speed U 10 ; steeper waves tend to occur under stronger winds.
In the lowest wind bin (blue) both C d and Z ch appear to increase with increasing H s /L p . The dependence of Z ch is roughly consistent with the previous parameterizations by Edson et al. (2013) and Taylor and Yelland (2001), based on observations in low to medium wind speeds. However, both C d and Z ch appear to remain flat or even decrease with increasing H s /L p in higher wind speeds (even within the same wind speed bin). Generally, the steepness H s /L p does not collapse the data very well and at least within our dataset. In the previous study (Hsu et al. 2019), which utilized the same ocean current observations, the wind stress was initially estimated from the vertical integral of two linear terms in the ocean current momentum equation, and was later corrected for the contribution from the neglected terms. In this study the MOM6-WW3 model was run to simulate the observed conditions using an assumed drag coefficient, and the drag coefficient was estimated for each float observation based on direct comparison between the observed and simulated ocean currents. This approach made it possible to utilize many data points behind the storm, to investigate the sea state dependence of C d in detail, and to investigate different surface wave impacts (Langmuir turbulence, air-sea momentum flux budget, wave-current interactions) on the upper-ocean currents and the C d estimates.
The estimated drag coefficient averaged over all TCs is consistent with the GFDL drag coefficient (Ginis et al. 2015), which is used as the assumed C d in the model simulations. Both our estimated C d and the GFDL C d are in the range of 2-3 3 10 23 for wind speeds of 25-55 m s 21 , and show weak wind speed dependence. This result is encouraging because the two drag coefficients have been obtained by completely different approaches: our C d is based on the upper-ocean current observations, while the GFDL C d is based on a large number of atmospheric model simulations of TCs. We find significant storm dependence of the estimated C d . We suspect that this is partly due to wind speed differences among different wind products under TC conditions. However, averaging over a large number of observations seems to minimize the impact of wind speed differences on the C d estimates.
One of the most significant findings in this study is that C d is significantly reduced by misaligned swell, that is, if the misalignment angle c between the dominant wave direction and the wind direction exceeds about 458. Since misaligned swell is common in the far front and in the front-left quadrant of the storm, C d tends to be lower in the same areas, and displays a distinct spatial distribution pattern as suggested by the previous study (Hsu et al. 2019). Our results also suggest that common occurrence of misaligned swell may be partially responsible for the overall reduction of C d under TC conditions. Misaligned swell also appears to turn the wind stress direction toward the swell direction from the wind direction.
In lower wind bins (25-30 m s 21 ) our estimates of the drag coefficient and the Charnock coefficient Z ch decrease with increasing wave age c p /U 10 . These trends are roughly consistent with previous observations made in low to medium wind speeds. However, in higher wind speed bins these trends disappear or even are reversed, that is, C d and Z ch may increase with increasing wave age c p /U 10 . We find that the input wave age c pi /U 10 and the wave steepness H s /L p do not distinguish different sea state effects at least within our dataset. These results are not conclusive because of the limited number of observational data. More observations are needed to refine such sea state-dependence analyses.
We find that the different surface wave effects (Langmuir turbulence, air-sea momentum flux budget, wave-current interactions) on upper-ocean currents/turbulence generally have minor impacts on the overall C d estimates. However, the wave-current interactions may have large impacts on the upper-ocean responses and the C d estimates in some conditions.
Finally, we have confirmed that our conclusions of the mean C d values as well as their dependence on sea states remain valid even if all the analyses are performed using a different wind product (see Supporting Information E in the supplemental material for more details).
Based on the findings discussed above, the two key conclusions from this study are as follows: 2) However, C d is more strongly dependent on sea states in the same wind speed range. In particular, this study clearly demonstrates that misaligned swell significantly reduces C d , a feature which is underestimated by current models of sea state-dependent drag coefficient.
The second conclusion strongly suggests that it is beneficial to implement an improved sea state-dependent parameterization of the drag coefficient in TC conditions. Some previous studies suggest that the drag coefficient parameterization should depend on storm relative locations. While our analyses do show spatial patterns of the swell misalignment angle c and other sea state parameters, it is unclear whether such an approach can fully capture significant sea state dependence. The simple quadrant analysis (front left, front right, etc.) is clearly insufficient.
Proposed alternative approaches are to estimate the sea state parameters based on the storm parameters (e.g., radius of maximum wind, maximum wind speed, storm translation speed) without running a wave model, such as using the effective wind duration parameter proposed by Hsu et al. (2019). We have found that this particular parameter is not able to distinguish the swell misalignment angle c and the resulting C d variation very well. Nevertheless, it may be possible to develop an alternative approach focused on accurately predicted c.
However, the sea state dependence of C d is expected to be complex. For example, the effect of misaligned swell may not depend on c alone, but may also depend on steepness and phase speed (relative to wind) of the swell. Therefore, it is likely more beneficial to couple a wave model and a hurricane (either atmosphere only or atmosphere plus ocean) model and introduce a sea state-dependent parameterization of C d based on the full wave spectrum.
Finally, developing such a parameterization requires better understanding of how misaligned swell reduces C d , based on more observational and theoretical/numerical studies, because existing models clearly underestimate this impact.
Brought to you by UNIVERSITY OF RHODE ISLAND | Unauthenticated | Downloaded 05/24/23 04:48 PM UTC
JULY 2022 Brought to you by UNIVERSITY OF RHODE ISLAND | Unauthenticated | Downloaded 05/24/23 04:48 PM UTC