A Novel Sea Surface Roughness Parameterization Based on Wave State and Sea Foam

: A wave state related sea surface roughness parameterization scheme that takes into account the impact of sea foam is proposed in this study. Using eight observational datasets, the performances of two most widely used wave state related parameterizations are examined under various wave conditions. Based on the different performances of two wave state related parameterizations under different wave state, and by introducing the effect of sea foam, a new sea surface roughness parameterization suitable for low to extreme wind conditions is proposed. The behaviors of drag coefﬁcient predicted by the proposed parameterization match the ﬁeld and laboratory measurements well. It is shown that the drag coefﬁcient increases with the increasing wind speed under low and moderate wind speed conditions, and then decreases with increasing wind speed, due to the effect of sea foam under high wind speed conditions. The maximum values of the drag coefﬁcient are reached when the 10 m wind speeds are in the range of 30–35m/s.


Introduction
The momentum transfer between the atmosphere and the ocean plays an important role in the evolution of weather and climate [1][2][3]. Parameterization of the momentum transfer across the air-sea interface is essential to the modeling of many air-sea interaction activities, such as tropical cyclones and ocean waves [4]. In the current applications, the air-sea momentum flux τ is usually estimated from the drag coefficient C d as follows: where ρ is the air density, u * is the friction velocity, and U 10 is the wind speed at 10 m elevation above the sea surface. The logarithmic wind profile law can be expressed as [5][6][7]: where κ = 0.4 is the von Kármán constant, and z 0 is the sea surface aerodynamic roughness, Ψ m is the stratification correction for the logarithmic wind profile, which is a function of the Obukhov length L, and the function of Ψ m can be found in Paulson [8] for unstable stratification and in Grachev et al. [9] for stable stratification, respectively. By combining Equations (1) and (2), the relationship between C d and z 0 is given as: Thus, there is an one-to-one correspondence between C d and z 0 under a certain stratification, specifying that z 0 specifies C d and vice versa. The sea surface aerodynamic roughness z 0 is widely used in the parameterization of the sea surface wind stress.
In current numerical models, C d and z 0 are often parameterized as the function of wind speed U 10 . In low and moderate wind conditions (U 10 ≤ 20 m/s), the results of many experiments show that C d increases linearly with wind speed [10][11][12][13]. Thus, the function of C d in low and moderate wind speed conditions can be expressed as [14]: By fitting the coefficients a and b to observational data, different results were obtained from different studies (Table 1); the functions of C d in low and moderate wind conditions from different research are qualitatively consistent, but differ significantly in values.

References a b
Kondo [15] 1.20 0.025 Smith and Banke [16] 0.63 0.066 Garratt [17] 0.75 0.067 Wu [18] 0.80 0.065 Large and Pond [19] 0.49 0.065 Donelan [20] 0.96 0.041 Geernaert et al. [12] 0.58 0.085 Yelland and Taylor [21] 0.60 0.070 Vickers and Mahrt [22] 0.75 0.067 Drennan et al. [23] 0.60 0.070 Guan and Xie [14] 0.78 0.065 Toffoli et al. [24] 0.96 0.060 Due to the lack of observational data in high wind speeds, the linear relationship between C d and U 10 in low and moderate winds has been extrapolated to high wind conditions in early applications, such as the modeling of tropical cyclones [25] and waves [26]. However, some recent experiments from both field and laboratory showed that C d tends to saturate [27,28] or decrease [29,30] with wind speed at extremely high wind speeds. Therefore, in many recent applications of tropical cyclone [31,32] and storm surge modeling [33], the increasing value of C d has been replaced by a constant that does not change with wind speed, or a value that decreases with increasing wind speed.
Several mechanisms of C d saturation at high wind speeds from different aspects have been proposed, and a summary of them can be found in Bryant and Akbar [34]. Many researchers ascribed the reduction or saturation of the C d to interface slipping and flattening accompanied by intense wave breaking at high wind speeds, which makes the wave steepness decrease or no longer increase, thereby affecting the aerodynamic roughness [35][36][37]. While some other researchers focused on the effect of sea foam on the momentum transfer process [38][39][40], the sea surface is covered by sea foam under high wind speed conditions, which changes the dynamics and thermodynamics of the air-sea interface. In addition to these two mechanisms, several other researchers explain the sea surface drag saturation from the unique airflow caused by breaking waves [41,42].
As the dependence of C d on wind speed varies significantly (Table 1), the drag coefficient might depend not only on the wind speed [43]. Based on the above mentioned mechanisms of C d saturation at high wind speeds, the dynamics and thermodynamics properties of the air-sea interface are crucial for the momentum transfer. Hence, it is convincible to parameterize the drag coefficient or the sea surface aerodynamic roughness through factors that describing the characteristic of the air-sea interface, i.e., wave age [44] and wave steepness [14].
Wave age and wave steepness are two of the most frequently used parameters to describe the air-sea interface and the development of wind wave. Wave age (β = c p /U 10 ) is defined as the ratio between spectral peak phase velocity c p and wind speed U 10 , or replace U 10 with friction velocity u * (β * = c p /u * ). Wave age β denotes the relative speed of wave to wind, the smaller the β, the lower the wave relative to the wind, and thus the more momentum transferred from the air to the sea. Wave steepness (δ = H s /L p ) is defined as the ratio between significant wave height H s and the wavelength at the spectral peak L p , δ denotes the physical roughness of the sea surface. In general, β describes the relative magnitude of wave speed and wind speed, while δ describes the characteristic of roughness.
Due to the importance of wave state on the momentum transfer across the air-sea interface, many wave parameter based schemes have been proposed to improve the parameterization of the momentum transfer [12,[45][46][47]. The dimensionless roughness z 0 /H s is often applied in the wave state related parameterization of the momentum transfer, as it has a stronger correlation with β and δ than the original C d and z 0 [48]. Smith et al. [11], Donelan et al. [46], and Drennan et al. [49] have proposed their function of z 0 /H s based on β or β * , respectively: These studies demonstrated a decreasing of the dimensionless roughness z 0 /H s with an increasing of wave age. On the other hand, Anctil and Donelan [50], Taylor and Yelland [51], and Takagaki et al. [28] have proposed their functions of z 0 /H s based on δ, respectively: z 0 /H s = 6.39 × 10 2 δ 6.76 , These studies demonstrate an increasing of the dimensionless roughness z 0 /H s with an increasing of wave steepness. The merits and limitations of both wave age based and wave steepness based sea surface roughness parameterization have been examined in several studies. Among them, the wave steepness based scheme proposed by Taylor and Yelland [51] (see Equation (9), hereafter TY01) and the wave age based scheme proposed by Drennan et al. [49] (see Equation (7), hereafter DN03) have received the most attention [48,52,53]. In general, the wave state related parameterizations present a better performance than the wind speed related bulk parameterizations, wave age based and wave steepness based schemes showed advantages in different wind or wave conditions, but none of them showed a good performance in all situations.
In addition to the wave state, sea foam also has a significant effect on the dynamics and thermodynamics properties of the air-sea interface. Under high wind conditions, the impact of sea foam on momentum transport cannot be ignored [54]. Owing to the lack of observational wave data under high wind conditions (U 10 ≥ 25 m/s), the aforementioned wave state related parameterizations have been proposed based only on the observational wave data under low to moderate wind speeds (U 10 ≤ 20 m/s). Note that, since the impact of sea foam on sea surface is minimal at low to moderate wind speeds [38,55], the effect of sea foam has not been included in these parameterizations implicitly.
In this study, we have evaluated the performance of two most widely used wave state related parameterizations (TY01 and DN03), using a combination of eight datasets including various wind and wave conditions. Based on the advantages and limitations of two schemes in different conditions, we propose a new wave state related parameteration scheme, by adding the effect of sea foam to the momentum transfer for existing schemes, which is verified to be suitable for low to extreme wind conditions (U 10 > 40 m/s).
The paper is organized as follows: Section 2 describes the observational datasets used to evaluate the performance of two wave state related parameterizations. Based on the different performances of two wave state related parameterizations under different wave states, a combination of them is proposed in Section 3. Section 4 introduces the effect of sea foam into the scheme presented in Section 3; thus, the new parameterization of sea surface roughness based on the wave state and sea foam is proposed. C d predicted by the new parameterization under high wind speed conditions is verified by the observational data in Section 5. Finally, Section 6 gives a summary of this study.

Datasets
To examine the performance of two most widely used wave state related parameterization: wave age based DN03 and wave steepness based TY01, eight observational datasets (published in tabular form) were used in this study. Wind stress in seven datasets was calculated using the direct eddy-correlation (EC) method [56] and the other dataset adopted the inertial dissipation (ID) method [57]. These datasets are described below, and a summary of them is given in Table 2.
a. Lake Ontario The Lake Ontario dataset was collected from the air-sea interaction experiment conducted in the western basin of Lake Ontario in the autumn of 1994 and 1995. A sonic anemometer was deployed on a 7.8 m-height bow mast to measure the wind fluctuations, which were used to calculate wind stress, and the sampling time of each run was 80 min (by pooling four consecutive 20-min averages groups to reduce the sampling error). Wave information was measured using a wave staff array. Here, we use the Lake Ontario data published by Anctil and Donelan [50].

b. AUSWEX
The Australian Shallow Water Experiment (AUSWEX) took place in the eastern basin of Lake George in 1997-2000 [58]. Two anemometer masts, accommodating wind probes, were mounted at 10-m height. Wind stress was calculated using the 21-Hz velocity data measured from an ultrasonic anemometer. Wave data were measured using eight wave probes. Here, we use the AUSWEX data published by Babanin et al. [59].

c. ERS Validation
The wind stress and wave data in this dataset were the validation data for the Grand Banks Earth Remote Sensing Satelite (ERS-1) Synthetic Aperture Radar (SAR) Wave Validation Experiment, which was collected from the scientific ship Hudson in the open North Atlantic. Wind data were measured using an anemometer system deployed on the bow of the ship, and the height of the system was 14 m. Wave data were measured using three wave buoys. Data used in this study were published by Dobson et al. [60].

d. SWADE
The data presented by Drennan et al. [61] were taken as part of the Surface Waves Dynamics Experiment (SWADE), which was conducted in 1990-1991 off the coast of Virginia. A 20-m swath ship was deployed to provide a high-resolution measurements near the air-sea interface [62]. Wind fluctuations were measured from an 12-m height anemometer, from which the wind stress was calculated, the sampling time was 17 min. Wave information was obtained using a wave staff array.

e. FPN
The North Sea Platform (FPN) experiment in 1985 was carried out on a platform located 65 km southwest of West Strand. Wind fluctuations were measured using a 33-m height sonic anemometer to calculate wind stress, and the sampling time was 30 min. Wave data were collected by a rider buoy located 800 m southwest of the platform, and were recorded on the platform. Data used in our study were released by Geernaert et al. [12] in tabular form.

f. HEXOS
The Humidity Exchange over the Sea (HEXOS) experiment was carried out on the Dutch research platform Meetpost Noordwijk (MPN) in the autumn of 1986. Wind fluctuations were obtained using a sonic and a pressure anemometer concurrently to calculated wind stress, height of them was 6 m, data collected from the pressure anemometer were adopted in this study, the sampling time for each run was 20 min. Wave data were collected by a rider buoy which was 150 m away from the platform. Here, we use the HEXOS data published by Janssen et al. [63].

g. RASEX
The Risø Air-Sea Exchange (RASEX) field experiment was performed at a shallowwater site near Denmark. In this experiment, wind fluctuation data were obtained from a 3 m height sonic anemometer, accompanied by the mean wind speed data collected from a cup anemometer located at 7 m, the sampling time for each run was 30 min. Wave data were gathered from the wave gauge near the tower. Data used here were obtained from Johnson et al. [64].

h. GOTEX
The Gulf of Tehuantepec Experiment (GOTEX) was carried out in February 2004. Data used in this study were presented by Romero and Melville [65], which were obtained from the National Science Foundation/National Center for Atmospheric Research (NSF/NCAR) C-130 aircraft. Vector winds were measured by the airborne detector at 25 Hz frequency, from which wind stress was calculated [66]. Frictional velocity u * was estimated from the lowest-height runs (about 40 m above the water surface) with a time average of 50 s. The sea surface elevation data were measured using a lidar system. In several datasets, the wavelength at the spectral peak L p was not measured directly. We calculate it using the dispersion relationship: where ω is the angular frequency, g is the gravitational acceleration, k is the wavenumber, and h denotes the depth of water. For deep water (h > L 2 , where L is the wavelength), L p can be calculated from: where T p denotes the period of wave at the spectral peak. If deep water conditions are not met, by substituting ω = 2π/T p and k = 2π/L p into Equation (11), L p can be calculated from T p and h. When both T p and f p (frequency of wave at the spectral peak) were not presented by the dataset, T p can be determined using the equations developed by Carter [67], which were derived from the Joint North Sea Wave Project (JONSWAP). For fetch limited seas: where X is the fetch in kilometers. For duration limited seas: where D is the duration in hours.
In addition to eight wind and wave datasets measured in low and moderate wind conditions, four datasets of C d in high wind speed conditions are also used in the validation of our new sea surface roughness parameterization in Section 5, two of them are field observations: Powell et al. [29] and Jarosz et al. [68], and the other two are laboratory observations: Donelan et al. [27] and Takagaki et al. [28]. Here, we make a brief introduction to them.
Powell et al. [29] measured the wind profile in tropical cyclone boundary layer using Global Positioning System, from the intercept and slope of the wind profile, C d and z 0 for winds up to 50 m/s are measured.
Due to the difficulties of direct stress measurements at high wind speeds caused by the spray droplets and the damages of winds to the instruments, Jarosz et al. [68] estimated the air-sea momentum transfer from the ocean side, namely the bottom-up method [27]. Using currents' observations recorded by the Acoustic Doppler current profiler during Hurricane Ivan, C d is calculated from: where ρ w and ρ are the density for water and air, respectively; f is the Coriolis parameter; U w and V w are the depth-integrated along and across the continental shelf current velocity components, respectively; U 10x is the along-shelf component of 10 m wind speed; and r is a constant resistance coefficient at the sea floor, which describes the degree of the bottom friction, and it usually ranges from 0.0001 cm/s to 0.1 cm/s. Using the bottom-up method, C d is estimated under different r for winds between 20 and 48 m/s. Donelan et al. [27] measured C d in laboratory conditions for winds up to 53 m/s using the Air-Sea Interaction Facility at the University of Miami, three methods were compared in the calculation of C d : momentum budget (MB), profile method (PM), and Reynolds stress (RS), the results from which were only slightly different. Tools for measuring stress include hot-film anemometry, digital particle image velocimetry (DPIV), and laser/line scan cameras for measuring the water surface elevation.
Using a high-speed wind-wave tank, Takagaki et al. [28] measured C d and z 0 for winds up to 64 m/s from wind velocity components collected by laser Doppler and phase Doppler anemometers; the eddy correlation method was utilized in their measurements to calculate C d and z 0 .

Evaluation of Two Wave State Related Parameterizations
The dimensionless roughness z 0 /H s of data points from eight datasets are plotted in Figure 1a against wave steepness δ, the curve of TY01 is also shown as the solid line.
It is shown that TY01 is able to describe the positive correlation between z 0 /H s and δ in general, but the data points are quite scattered.
For comparison, we plot the same data points using the wave age scaling in Figure 1b, i.e., z 0 /H s versus β * . The curve of DN03 provides a better prediction of the dimensionless roughness z 0 /H s than TY01, the data points are more concentrated near the curve than in Figure 1a. Although the plots of dimensionless roughness z 0 /H s for all data points show the overall performance of two parameterizations, it is more instructive to test how they predict the drag coefficient. z 0 can be converted to C d by Equation (3). A comparison between measured and predicted C d has been made for each dataset, and the results of TY01 and DN03 are presented in Figures 2 and 3, respectively. Note that the data points that fall within the 90% confidence regions are denoted as black points. The 90% confidence regions for datasets using the EC method are calculated based on the sampling errors ε [69], where the sampling errors of six EC datasets can be calculated following Donelan [70]: where Υ is the sampling time (s), U is the mean wind speed for an experiment, and z is the height of the anemometer above the water level. The sampling errors of eight datasets are summarized in Table 3. It is worth mentioning that the wind stress in the ERS Validation dataset was calculated using the ID method, Equation (16) is not applicable, and, following Drennan et al. [52], we assume an error equal to the mean sampling error of the EC data (25.77%). Similarly, data from GOTEX were collected from the aircraft, and the measuring instrument and post-processing method were inconsistent from other datasets. Equation (16) is suitable mainly for traditional platforms, i.e., buoy and tower. Thus, the sampling error for GOTEX dataset was also assumed as the mean sampling error for the EC data (25.77%). In Figures 2 and 3, the 90% confidence regions are shown as the areas between the dotted lines, and the slope of the upper and the lower boundary line is 1 + ε and 1/(1 + ε), respectively. To evaluate the performance of TY01 and DN03 quantitatively, P 90 was defined as the percentage of data points that fall within the 90% confidence regions. The normalized bias (NB) is defined as: and the normalized root-mean-square-error (NRMSE) is defined as: where X obs is the observation, and X mod is the corresponding value calculated from parameterization schemes [71]. In addition, P 90 , NB, and NRMSE predicted by TY01 and DN03 for each dataset are shown in Table 4, also shown are the mean β, mean β * , and mean δ for each dataset. From Table 4, we can see that the correlation between P 90 and NB or NRMSE is strong, and datasets with larger P 90 tend to have smaller NB and NRMSE, and datasets in which TY01 performs better under P 90 are consistent with that under NRMSE. Consid-ering that P 90 is consistent with NB and NRMSE qualitatively, and has the advantage of being able to take into account the sampling error of each dataset, we mainly focus on P 90 in the following analysis. The sampling errors of the ERS validation dataset and the GOTEX dataset were assumed as the mean sampling error of six other datasets. We first consider the results predicted by TY01 shown in Figure 2 and Table 4. TY01 is seen to work well for the AUSWEX, HEXOS, RASEX, and GOTEX datasets with a P 90 larger than 0.65, but C d measured in ERA Validation, SWADE, and FPN datasets was poorly predicted with a P 90 less than 0.4, especially in the FPN dataset (P 90 = 0.1638). As we can see from Figure 2, C d in ERS Validation and FPN datasets were extremely overpredicted by TY01, it is worth noticing that the mean β * of ERS Validation and FPN datasets were the largest two among eight datasets (both larger than 20), corresponding to a mature wave field. Moreover, TY01 underpredicted C d from AUSWEX and RASEX datasets, whose mean β * was the smallest two among eight datasets. The performance of TY01 shows an obvious sensitivity to β * ; for datasets having a larger β * , TY01 tend to overpredict C d from them; but, for datasets having a smaller β * , C d from them was underpredicted.
The results of DN03 were shown in Figure 3 and Table 4. The overall performance of DN03 is better than TY01. The results of DN03 from ERS Validation, SWADE, and FPN datasets are much better than TY01, but C d measured in AUSWEX, HEXOS, and RASEX was poorly predicted and worse than TY01. For datasets in which DN03 performs well, the mean β * was seen to be large (20.89, 18.88, 27.43, and 17.69 for ERS Validation, SWADE, FPN, and GOTEX, respectively), and in two datasets that have a smaller β * (7.54 for AUSWEX and 12.68 for RASEX), the performance of DN03 is quite worse. Therefore, the performance of DN03 also shows a sensitivity to β * .
In order to analyze the applicability of TY01 and DN03 in different conditions, we examine the sensitivity of their performance to β, β * , and δ. Here, we use TY01_in to denote the data points predicted by TY01 that fall within the 90% confidence regions, corresponding to those data accurately predicted by TY01; and TY01_out to denote the data points predicted by TY01 that fall outside the 90% confidence regions, corresponding to those data that are not accurately predicted by TY01. DN03_in and DN03_out are the same, but for data points predicted by DN03. Table 5 shows the mean β, mean β * , and mean δ of TY01_in, TY01_out, DN03_in, and DN03_out.
The mean β and β * of TY01_in is much smaller than that of TY01_out, demonstrating that TY01 tends to have better performance at younger wave conditions. The mean δ of TY01_in is close to the mean δ of TY01_out, indicating that the performance of TY01 is not sensitive to δ. The difference of the mean β between DN03_in and DN03_out is not as obvious as between TY01_in and TY01_out, but the difference of the mean β * between DN03_in and DN03_out is non-negligible. The difference of the mean δ between DN03_in and DN03_out is not obvious, demonstrating that the performance of DN03 is also not sensitive to the wave steepness. Considering that the performance of TY01 and DN03 is both sensitive to β * , to further investigate the sensitivities of the performance of TY01 and DN03 to β * , we divide the 471 data from eight datasets into 10 groups of roughly equal numbers (47 or 48 per group) according to β * from low to high, and calculate the P 90 of each group, the results are shown in Table 6. Changes in performance of TY01 and DN03 with β * are clearly demonstrated, when β * exceeds 16, the performance of TY01 drops significantly; when β * is smaller than 10, the performance of DN03 is relatively poor. Considering the different performance of TY01 and DN03 in different conditions, it is reasonable to combine them by using TY01 in small β * conditions and using DN03 in large β * conditions. Another issue is the choice of the demarcation point between TY01 and DN03, since the datasets used in this study do not cover all wind and wave conditions, and there are inconsistencies between datasets due to different observation and processing methods, we cannot determine the demarcation points arbitrarily as the point where the performance of DN03 exceeds TY01. Therefore, we use the δ − β * relationship derived from Toba's [72] 3/2 power law to determine the demarcation points between TY01 and DN03. The well-known 3/2 power law is given as: where H * = gH s /u 2 * and T * = gT s /u * are non-dimensional significant wave height and period, and B = 0.062 is a constant. The 3/2 power law has been verified by many studies [64,73,74], which is suitable for low to extreme wind conditions [48]. Multiplying Equation (19) by 2πu 2 * /g 2 T 2 p , we get: by using the relation between significant wave period T s and peak wave period T p [75,76]: and by calling the relation c p = gT p /2π, Equation (20) can be rewritten as: By combining Equation (7) (function of DN03), Equation (9) (function of TY01), and Equation (22), we work out that the curves of TY01 and DN03 intersect at β * = 15.21. According to the above inference, β * = 15.21 is selected as the demarcation point between TY01 and DN03, TY01 is adopted when β * < 15.21, and DN03 is adopted when β * ≥ 15.21: To verify the validity of the combination of TY01 and DN03 given in Equation (23), Figure 4 plotted a comparison between the measured C d and the corresponding values predicted by Equation (23) as in Figures 2 and 3. By comparing Figures 2 and 4, we can see that the performance of the combined scheme is much better than that of TY01, especially in ERS Validation and FPN datasets. By comparing Figures 3 and 4, the improvement of the combined scheme compared to DN03 mainly comes from the RASEX dataset; most of the RASEX data overestimated by DN03 have been improved in the combined scheme. We further compared the P 90 , NB, and NRMSE predicted by TY01, DN03, and the combined scheme for the total eight datasets (Table 7); the results show that the performance of the combined scheme is much better than TY01 in P 90 and NRMSE, and slightly better than that of DN03, NB predicted by the combined scheme is slightly worse than DN03. Considering that NB mainly describes the overestimation or underestimation of the prediction, and can be offset if both overestimation and underestimation exist, while P 90 and NRMSE are the key parameters to show the overall performance; the results in Table 7 prove that the performance of the combined scheme is better than TY01 and DN03.

Effect of Sea Foam
TY01 was developed using three datasets: HEXOS, RASEX, and Lake Ontario, and DN03 was developed using the pure wind sea subsets of five datasets: AGILE (measured from the 15-m research vessel AGILE) [77], FETCH (Flux, sea state and remote sensing in conditions of variable fetch) [78], HEXOS, SWADE, and WAVES (Water-Air Vertical Exchange Study) [79]; these datasets were collected under low and moderate wind conditions (U 10 ≤ 20 m/s). Compared to low and moderate wind conditions, a significant change in high wind conditions (U 10 ≥ 25 m/s) is the generation of sea foam due to intense wave breaking, which plays an important role in the leveling off or decrease of C d and z 0 . Since the effect of sea foam on sea surface roughness is minimal at low to moderate wind speeds [38,55], the effect of sea foam was not implicitly included in the proposing of TY01 and DN03, and an introduction of the effect of sea foam to TY01 and DN03 will enhance their applicability for high wind speed conditions.
A semi-empirical model is proposed by [55] to estimate the influence of sea foam on aerodynamic roughness. Their model treats the effective air-sea aerodynamic roughness (z e f f ) as the weighted sum of two parts: one is the foam-free (z n ) part and the other is the foam-covered (z f ) part. The average z e f f under area S is assumed as follows: Here, S = S n + S f is the total area, in which S n and S f are the foam-free and foam-cover areas, respectively. Thus, Equation (24) can be rewritten as: by defining α f = S f /S as the fractional foam coverage, we obtain: The fractional foam coverage α f is highly related to U 10 [80]. The function between α f and U 10 can be approximated from the observational data as in Holthuijsen et al. [80]: with α = 0.00255, ζ = 0.166, and γ = 0.98. To demonstrate the different patterns of α f in different situations, a universal dimensionless form of Equation (27) is given as: 10 is the saturation speed, defined as the value where the difference between α f and its saturation limit α f = 1 is less than 2%. The curve of foam coverage α f versus U 10 from Equation (28)  , z e f f [29], it is assumed that the minimum value for C d = 0.0017 or, z e f f = 0.0003 m is reached at the same wind speed U 10 = 48 m/s (see Figures 2 and 3 in Golbraikh and Shtemler [55]). Because the relation between C d and U 10 in laboratory conditions is quite different from that of the open ocean, Golbraikh and Shtemler [55] suggested a different minimum value for z e f f in laboratory conditions, which is z e f f = 0.0028 m. Then, we adopt U 10 = 48 m/s as the saturation velocity, and the minimum value of z e f f = 0.0003 m as the foam-covered aerodynamic roughness z f in Equation (26) for open ocean conditions, and the minimum value of z e f f = 0.0028 m for laboratory conditions. As the effect of sea foam was not implicitly included in the proposing of TY01 and DN03, the aerodynamic roughness predicted by Equation (23) can be taken as the foam-free aerodynamic roughness z n in Equation (26), substituting Equation (23) into Equation (26), a new parameterization of sea surface roughness including the impact of sea foam is obtained: where z 0 is the aerodynamic roughness, δ is wave steepness, β * is wave age, H s is the significant wave height, α f is the foam coverage (calculated from Equation (28)), and z f is the foam-covered aerodynamic roughness (taken as 0.0003 m for open ocean conditions, and 0.0028 m for laboratory conditions in this study). By combining TY01 and DN03 in the form of a piecewise function, the new proposed parameterization is able to make better predictions of z 0 in various wind and wave conditions. By adding the impact of sea foam, the predictions of z 0 in high wind speed conditions are improved.

Validation and Discussion
As aforementioned, the proposed parameterization calls TY01 and DN03 according to different wave ages; its comparison against observations in low and moderate wind conditions have been made in Section 3. In this section, we will make a brief validation on the behavior of C d predicted by the proposed parameterization under high wind speed conditions. Specifically, datasets from several recent experiments [27][28][29]68] with observations under high wind speed conditions are compared with the new parameterization.

Estimation of H s
These observational data were presented in the form of C d vs. U 10 , in order to compare the proposed parameterization with these observational data; it is essential to parameterize H s with U 10 .
Several schemes were proposed for the parameterization of H s . From the formulas for fully developed wave field in deep water, Taylor and Yelland [51] proposed a parameterization of H s : According to Equation (30) These schemes all reveal the monotonically increasing of H s with U 10 , and this trend has been verified in low and moderate wind speeds, their applicability in high wind speed conditions is doubtful. The plots of H s versus U 10 of the three schemes above are shown in Figure 6a, the values of H s from three schemes are relatively reasonable at low and moderate wind speeds; however, as the wind speed increases, H s becomes unreasonably large, the values of H s calculated from three schemes all exceed 50 m at U 10 = 60 m/s, which are obviously unreasonable. However, accurate prediction of H s under high wind speeds requires the help of numerical models. Considering that our purpose is only to get the brief relationship between H s and U 10 , we simply add a threshold of 21 m to H s to replace the unreasonably large value under high wind speeds (Figure 6b). The value of 21 m comes from the largest H s measured by the radar altimeter onboard the Jason 2 satellite (http://cersat.ifremer.fr/user-community/news/item/346-record-breaking-waveheights-and-periods-in-the-north-atlantic, accessed on 11 February 2021) which is 20.1 m; here, we round it to 21 m.

Validation of the Proposed Parameterization
In order to show how sea foam affects our results, the comparison between the curves of our parameterization without the effect of sea foam (Equation (23)) and the field observations from Powell et al. [29] and Jarosz et al. [68] is presented in Figure 7. Figure 7a-c denotes the different relations from H s estimated from Taylor and Yelland [51], Fairall et al. [81], and Wang et al. [82], respectively. For the curves of β * < 15.21, δ has been converted to β * using the δ − β * relationship derived from Toba's [72] 3/2 power law (Equation (22)). From Figure 7, we can see that the curves of C d from Equation (23) can not reproduce the decreasing of C d at high winds. The effect of sea foam can be seen from the comparison between Figures 7 and 8.
Given that z f in Equation (29) is taken as different values for field and laboratory conditions, we compared the new parameterization with field and laboratory observations separately. Figure 8 shows the comparison between C d predicted by the new parameterization (Equation (29)) under different wave ages and the field observations from Powell et al. [29] and Jarosz et al. [68]. From Figure 8, we can see that C d predicted by the new proposed parameterization using different H s schemes are generally consistent. C d increases with wind speed in the range of 0-30 m/s, the maximum values are reached at about 30∼35 m/s, then decreases at the wind speed about 35∼45 m/s under the effect of sea foam, for wind speed larger than 45 m/s, the values of C d do not change much. By comparing Figures 7 and 8, the effect of sea foam is obvious, by adding the sea foam, our parameterization can reproduce the reduction of C d at U 10 > 30 m/s, which is closer to the observations.
Results from H s schemes proposed by Taylor and Yelland [51] and Fairall et al. [81] (Figure 8a,b, respectively) do not show much differences, but the results from Wang et al. [82] (Figure 8c) are different from the other at low and moderate wind speeds, in which the values of C d are larger than the other two, especially for the younger wave. The difference is caused by the intercept of the formula proposed by Wang et al. [82] (see Equation (32)), when U 10 is close to zero, H s still has an initial value, given that young wave fields generally do not correspond to low wind speeds; this difference is not obvious in practice.  (23)) under different wave ages using H s estimated from (a) Taylor and Yelland [51]; (b) Fairall et al. [81]; and (c) Wang et al. [82] and the field observations.
The curves of the new proposed parameterization shown in Figure 8 can cover the range of the field observational data well, and the scatter of the observations can be explained as the effect of wave state. The reduction of C d under high wind speeds is successfully reproduced by the new proposed parameterization, C d predicted by the new proposed parameterization reach the maximum values in the wind range of 30∼35 m/s, which is consistent with the field measurements in Jarosz et al. [68] and Powell et al. [29]. The maximum value in Jarosz et al. [68] with the resistance coefficient of 0.1 cm/s (∼3.7 × 10 −3 ) is close to the maximum value of β * = 9 in Figure 8a,b, and is between β * = 9 and β * = 6 in Figure 8c. Furthermore, compared with the curves in Jarosz et al. [68] (cf Figures 2 and 3 therein), our parameterization provides C d values for U 10 > 50 m/s, while curves in Jarosz et al. [68] did not, considering that conditions with U 10 larger than 50 m/s are common in tropical cyclones, our parameterization is suitable for the usage in tropical cyclone modeling and storm surge modeling. Since the simultaneous wave state was not measured by Jarosz et al. [68], we cannot compare the predictions of C d with the observations directly. Figure 9 shows the comparison between C d predicted by the new parameterization under different wave ages and the laboratory observations from Donelan et al. [27] and Takagaki et al. [28]. The laboratory measurements do not show a decreasing trend under high wind speeds, their C d tend to saturate at wind speeds larger than 35 m/s. The difference between the field and laboratory measurements can be expected due to significant differences in fetch [35]; in addition, in hurricane conditions, the wave field is dominated by swell generated in the high wind areas, but it will not be reproduced under laboratory conditions [83]. Consistent with the laboratory measurements, C d predicted by the new parameterization also shows a saturation at U 10 > 40 m/s, the saturation values of the observations match the predicted C d well, both of them are very close to C d = 0.0024. For U 10 < 30 m/s, the values of observations concentrate near the curve of a larger wave age; considering that wave age is negative related to wind speed, and lower wind speed usually corresponds to a larger wave age, this result is reasonable. Observations from Donelan et al. [27] are slightly lower than that predicted by our parameterization, especially for MB and PM methods, this slightly difference is caused by the calculation method and the measuring instrument, i.e., the RS method uses the stress data directly measured from an x-film anemometer, the PM method uses wind speed data measured from the hot-film anemometry, and the MB method uses the bottom stress from DPIV and surface elevation from laser/line scan cameras to calculate C d .
Although the proposed parameterization can reasonably explain the behavior of the observational data, it should be pointed out that the values predicted by the new parameterization have not been compared with the observational data directly due to the lack of simultaneous wave state measurements under high wind speed conditions. Thus, more field and laboratory experiments containing simultaneous wind and wave state measurements are needed to further verify the performance of the new parameterization, and to investigate the mechanism of momentum transfer across the air-sea interface.  β* = 6 β* = 9 β* = 12 β* = 15 β* = 18 β* = 21 β* = 24 β* = 27 β* = 30 Figure 9. Comparison between C d predicted by the new parameterization under different wave ages using H s estimated from (a) Taylor and Yelland [51]; (b) Fairall et al. [81]; and (c) Wang et al. [82] and the laboratory observations.

Comparison with Other Parameterizations
In this section, the performance of our parameterization has been compared with three different parameterizations, these parameterizations have been proposed for the calculation of C d in high wind speed conditions, and the saturation of C d has been dealt with different method.
Based on the work of Powell [84] and Garratt [17], Luettich and Westerink [85] offered a formulation that divides the tropical cyclone into three sectors and calculated C d accordingly; this formula has been used in the ADvanced CIRCulation (ADCIRC) storm surge model (here, we denote it as ADCIRC). For the right sector of a storm: for the rear sector of a storm: and for the left front sector of a storm: Using more than 6000 near-surface flux measurements collected from low-flying aircrafts, Andreas proposed a parameterization for low-to-high winds (here, we denote it as A12): C d can be calculated from C d = ( u * U 10 ) 2 . According to the dependence of wind speed-C d relation on swell, Holthuijsen et al. [80] proposed a parameterization for different swell conditions (here, we denote it as H12): For no swell, opposing swell, and following swell, a = 1.05, b = 1.25, c = 1.4, d = 2.3, and e = 10; for cross swell a = 0.7, b = 1.1, c = 6, d = 8.2, and e = 2.5.
The comparison between the above-mentioned three parameterizations and our new parameterizations, along with two field observations, are presented in Figure 10 because the curves of our parameterization only show a little difference for different H s parameterization as shown in Figures 8 and 9, we only plot the curves based on the H s parameterization of Wang et al. [82]. The results in Figure 10 show that the saturation of C d are presented in different forms, all three forms of ADCIRC take C d as constants for U 10 > 45 m/s, but the values are different, the maximum value of C d predicted by ADCIRC left front is about 0.0057, which is much larger than that observed in field experiments; another problem for ADCIRC parameterization is that their values of C d are not continuous at their demarcation point of their formula, such as U 10 = 30 m/s, U 10 = 35 m/s, and U 10 = 45 m/s, which are unreasonable physically. The curve of A12 is smooth, but the reduction of C d at U 10 > 35 m/s has not been reproduced by their formula. C d predicted by H12 under cross swell conditions reaches a maximum value of about 0.0053 at U 10 ≈ 35 m/s, then de-creases rapidly, the predicted C d is smaller than 0 when U 10 > 54 m/s, this is also incorrect. In general, these three schemes have deficiencies in different aspects, our parameterization has presented the most reasonable results.

Conclusions
An accurate estimate of momentum transfer across the air-sea interface is vital for atmospheric, oceanic, and surface wave prediction models. Compared with parameterization of momentum flux based on wind speed, parameterization based on wave state can describe the nature of the air-sea interface more directly. Wave age (β = c p /U 10 , or β * = c p /u * ) and wave steepness (δ = H s /L p ) are two of the most frequently used parameters to describe the air-sea interface and the development of wind wave. Using eight observational datasets, the performances of two most widely used wave state related parameterizations: TY01 and DN03, are examined under various wave conditions. TY01 shows a better performance for the younger waves (smaller β * ), while DN03 is more suitable for wave fields with medium or large wave age. Hence, we use a combination of them to get a better performance under various wave conditions: for β * < 15.21, TY01 is adopted; and, for β * ≥ 15.21, DN03 is adopted. The demarcation point β * = 15.21 is selected from the δ − β * relationship derived from Toba's [72] 3/2 power law (see Equation (22)). Considering that TY01 and DN03 were developed using observational data under low and moderate wind speed conditions (U 10 ≤ 20 m/s), the effect of sea foam was not included explicitly or implicitly in the proposing of TY01 and DN03. By introducing the effect of sea foam into the scheme presented in Section 3 (see Equation (23)), a new parameterization of sea surface roughness based on the wave state and sea foam is proposed (see Equation (29)). C d predicted by the new parameterization increases with wind speed in the range of 0∼30 m/s; the maximum values are reached at about 30∼35 m/s and then decrease at the wind speed about 35∼45 m/s under the effect of sea foam; its behavior is also supported by the field observations from [29,68]. The saturation values of C d in laboratory measurements from [27,28] is also reproduced by the new parameterization.
Due to the vital role of wave state and sea foam on the momentum transfer across the air-sea interface, the new proposed sea surface roughness parameterization is suitable for the coupled atmosphere-ocean-wave modeling systems. Furthermore, as the effects of sea foam are included in the presented parameterization, it is also applicable for the modeling of some severe air-sea interaction activities accompanied with extreme winds, such as tropical cyclones, and the wave modeling of storm surge.
Finally, it should be emphasized that, due to the lack of simultaneous wave state measurements under high wind speed conditions, the values predicted by the new parameterization have not been compared with the observational data directly. Thus, more field and laboratory experiments containing simultaneous wind and wave state measurements, especially for high wind speed conditions, are needed to further verify the performance of the new parameterization, and to investigate the specific mechanism of air-sea interaction. However, although a direct comparison between the new parameterization with the observational data at high wind speeds is difficult, assessing it in the numerical weather prediction system is more realistic. It is our plan to implement the new parameterization in numerical models, including large-eddy simulations and coupled atmosphere-wave models, and to evaluate the performance of our parameterization from the model results.