Assessment of surface drag coefficient parametrizations based on observations and simulations using the Weather Research and Forecasting model

Abstract The drag coefficient is important in meteorological studies of the boundary layer because it describes the air–sea momentum flux. Eight drag coefficient schemes were assessed. These parametrizations were compared taking into account data from in situ and laboratory observations. The drag coefficients determined using three schemes were consistent with the level-off phenomenon, supported by the results of laboratory studies. The drag coefficient determined using one scheme decreased at wind speeds higher than approximately 30 m s−1, in agreement with indirect measurements under typhoon conditions. In contrast, the drag coefficients determined using the other four schemes increased with wind speed, even under high wind regimes. Sensitivity tests were performed using simulations of two super typhoons in the Weather Research and Forecasting model. While the typhoon tracks were negligibly sensitive to the parametrization used, the typhoon intensities (the maximum 10-m wind speed and the minimum sea level pressure), sizes, and structure, were very sensitive to it.


Introduction
A bulk formula is often used to represent the exchange of momentum flux (also known as wind stress) between the atmosphere and the ocean surface. The bulk formula is given by where τ is the air-sea momentum flux, ρ is the density of the air, U 10 is the 10-m wind speed, and C d is the drag coefficient. The extrapolation of C d derived from field measurements made in weak wind regimes (assuming that C d continuously increases with wind speed) has been shown not to be robust either in theory (Andreas 2004;Emanuel 2003), laboratory experiments (Alamaro 2001;Donelan et al. 2004;Troitskaya et al. 2012), or field studies (Bell, Montgomery, and Emanuel 2012;Black et al. 2007;Jarosz et al. 2007;Potter et al. 2015;Powell, Vickery, and Reinhold (1) = C d U 2 10 , 2003). It is believed that C d stops increasing for strong winds (generally with wind speeds greater than 30 m s −1 ). Various parametrizations of C d have been proposed in recent years, but the parametrizations differ greatly because they have been developed using a variety of methods and datasets. A comprehensive comparison of the parametrizations would therefore be useful for those wishing to use C d in practice. The air-sea momentum flux plays a crucial role in the intensification and maintenance of a tropical cyclone. Potential intensity theory and the results of numerical models presented by Emanuel (1986Emanuel ( , 1995 suggest that the square of the maximum wind speed (V 2 max ) and the minimum central pressure (P min ) of a tropical cyclone in a given environment will be dependent on the ratio of the enthalpy exchange coefficient (C k ) and C d , and that the C k /C d ratio in a real typhoon should not be less than 0.75. Other modelling studies (Bryan 2012(Bryan , 2013Bryan 1989), and the 'rapid radiative transfer model' longwave scheme (Mlawer et al. 1997) were implemented in both domains. The Kain-Fritsch cumulus scheme (Kain 2004) was used only in D01, and no cumulus parametrization was adopted for D02. A simple mixed-layer ocean model (Pollard, Rhines, and Thompson 1973) was coupled with the atmospheric model to provide updated SSTs.
The initial conditions and lateral boundary conditions in the WRF model were obtained from the NCEP 'FNL' (final) 'operational global analysis' data-set, with a 1° × 1° grid and at intervals of six hours. The wave-age-related parametrization of C d was examined using peak-wave-period data with a resolution of 0.5° (from the US Environmental Research Division's Data Access Program; http://oos.soest.hawaii. edu/erddap/griddap/NWW3_Global_Best.html) and providing these to the WRF model every six hours to allow the wave age to be calculated. The peak-wave-period data were taken from the WaveWatch III Global Wave Model, implemented by the University of Hawaii. The wave model was forced using the U.S. NCEP 'Global Forecast System' wind data, and it was initialized daily, making the peakwave-period data consistent with the initial and lateral boundary conditions of the WRF model.
The experiments were begun at 1800 UTC 4 November 2013 for Typhoon Haiyan and at 1800 UTC 28 March 2015 for Typhoon Maysak, and were integrated for 72 h, covering the development and rapid strengthening stages of the typhoons. The hourly results from D02 were used in the analysis described below.

Description of the selected parametrizations
In a neutrally stratified atmosphere, the surface C d can be expressed in terms of the momentum roughness length (z 0 ) as where the von Karman constant k is taken to be 0.4 and the reference height z is typically 10 m.
The default option for the z 0 parametrization in the WRF model is based on the Charnock relationship given by with the Charnock coefficient α = 0.0185 (Wu 1980), where g is acceleration due to gravity and u * is the friction velocity. The second term on the right-hand side of Equation (3) represents the smooth-flow component, supported by viscous shear (Edson et al. 2013).
(2) and Rotunno 2009;Green and Zhang 2013) have produced results that are consistent with the potential intensity theory to some degree, and it has been accepted that V max ~ (C k /C d ) 1/2 . Green and Zhang (2013) conducted five physical sensitivity experiments using three different momentum flux parametrizations and two different moist enthalpy flux parametrizations, using version 3.4.1 of the Weather Research and Forecasting (WRF) model. They found that the surface fluxes significantly affected a simulation of Hurricane Katrina (2005) and that in particular C d affected the pressure-wind relationship and the maximum surface wind radius of the hurricane. Recently, it has been shown that C d can be optimized using data assimilation (e.g. 4D Variational Data Assimilation), and some new models can be obtained through data assimilation Peng and Li 2015;Peng, Li, and Xie 2013).
In the study, eight previously published C d parametrizations were compared with each other and with observational data. The sensitivities of simulated typhoons to the different parametrizations were then examined using version 3.6 of the WRF model.
The configuration of the WRF model and the experimental design are presented in Section 2. The eight C d parametrizations are described in Section 3. A comparison and validation of the parametrizations with the observations, as well as analyses of the modelling results for two super typhoons, are provided in Section 4. And finally, the study is summarized in Section 5.

Method
The 'Advanced Research' core of the WRF model, version 3.6, was utilized to simulate typhoons Haiyan and Maysak using the selected C d parametrizations. A two-domain configuration with two-way nesting was used in the WRF model. The horizontal grid resolution in the outer domain (D01) was 18 km, and in the inner domain (D02) was 6 km. For the simulation of typhoon Haiyan (2013), D01 was from 0° to 20°N and 115°E to 155°E, and D02 was from 4°N to 14°N and 120°E to 150°E. For the simulation of typhoon Maysak (2015), D01 was from 0° to 20°N and 120°E to 170°E, and D02 was from 5°N to 15°N and 134°E to 156°E. The time steps used in D01 and D02 were 60 and 20 s, respectively. There were 44 levels in the vertical direction, and 12 of these levels were inhomogeneous and below 3 km. The model extended up to a height of 50 hPa (approximately 30 km). The microphysics scheme presented by Lin, Farley, and Orville (1983), the revised MM5 Monin-Obukhov surface layer scheme (Beljaars 1995;Dyer and Hicks 1970;Paulson 1970;Webb 1970;Zhang and Anthes 1982), the BouLac planetary boundary layer scheme (Bougeault and Lacarrere 1989), the Dudhia shortwave scheme (Dudhia Taylor and Yelland (2001) proposed that z 0 can be predicted from the significant wave height (h s ) and the wave steepness (h s /L p ) so that where L p is the peak wavelength. For a developed sea, h s and L p are determined from U 10 according to and Moon et al. (2007) analyzed 10 tropical cyclones simulated using a coupled wave-wind model, and then developed a new parametrization of z 0 in the form Davis et al. (2008) used the experimental results of Donelan et al. (2004) to derive an alternative formula for z 0 , where and The schemes presented by Taylor and Yelland (2001) and Davis et al. (2008) are both optional in the WRF model, but the z 0 has an upper limit of 2.85 × 10 −3 m. In this study, these two schemes without the upper limits, as well as the Taylor and Yelland scheme with the upper limit, were examined. Foreman and Emeis (2010) and Andreas, Mahrt, and Vickers (2012) suggested that the traditional expression u * = aU 10 , transformed from C d = (u * /U 10 ) 2 , requires an extra constant b, to give Andreas, Mahrt, and Vickers (2012) calculated factors a and b by analyzing almost 7000 air-sea momentum flux measurements, and they found a new quantitative relationship between u * and U 10 : We used the logarithmic wind profile to obtain the formation of z 0 : Edson et al. (2013) used four oceanic field experiment datasets to investigate the air-sea momentum flux, and modified the Charnock coefficient α to α = 0.0017U 10 −0.005 with an upper limit of 0.028, This parametrization has been added to the more recent version (3.7) of the WRF model. Liu, Guan, and Xie (2012) proposed a wave-related parametrization of z 0 , based on the theoretical results of the Scientific Committee on Oceanic Research workgroup 101 (Jones and Toba 2008), taking the effects of ocean waves and sea spray on α into consideration: where β * = c p /u * is the wave age, c p is the phase speed at the peak of the wave spectrum, calculated using c p = gT p /(2n) in the WRF model, and ω = min (1.0, 1.6/u * ) is a correction factor indicating the influence of sea spray on C d .

Comparison and validation of the selected parametrizations
The parametrizations of z 0 that were tested are referred to using the author and year of the publication in which the parametrizations were described. We thus use (7a) u * = aU 10 + b.
(7b) u * = 0.0283U 10 + 0.00513, U 10 ≤ 9 m s −1 0.0583U 10 − 0.243, U 10 > 9 m s −1 . (8) (10) = min[(0.0017U 10 − 0.005), 0.028]. (11) , rapidly (at almost the same rate) in the Charnock1955 and Davis2008 scheme. The value of C d in the TY2001 scheme increased at a moderate rate and the C d became smaller than in the Davis2008 scheme at wind speeds of 40 m s −1 and higher. The C d tended to level off and approach roughly 3 × 10 −3 at high wind speeds in both the Andreas2012 and the Moon2007 scheme. The upper limit of 2.85 × 10 −3 m on z 0 caused the C d in the TYlimit scheme to remain at approximately 2.4 × 10 −3 at wind speeds higher than 20-30 m s −1 . The wave-age-related Liu2012 scheme was distinctly different from the other schemes in that C d reached a maximum at wind speeds of 25-33 m s −1 and then decreased linearly as the wind speed increased further. A lower wave age in the Liu2012 scheme corresponded with a higher maximum C d and weaker surface winds. In the real simulation of Typhoon Haiyan, the wave age was found to decrease quickly as the wind speed increased, causing C d to increase quickly in weak and moderate winds and then slowly decrease in strong winds (figures not shown).
The C d values determined using the eight parametrizations and observational data are also compared in Figure 1. The C d values determined using the Andreas2012, Moon2007, and TYlimit schemes were very close to the C d values obtained from laboratory measurements made by Donelan et al. (2004) and Troitskaya et al. (2012), particularly at high wind speeds (30-52 m s −1 ). The C d only decreased at wind speeds higher than 25-33 m s −1 in the Liu2012 scheme, and this agreed with field observations reported by Powell, Vickery, and Reinhold (2003) and Jarosz et al. (2007). The Liu2012 scheme also partly attributed the scatter in the observed C d distribution at specific wind speeds to the different wave states during the measurements (Liu, Guan, and Xie 2012). However, the observations in Figure 1 were obtained in a limited range of wind speed (less than 52 m s −1 , roughly). Bell, Montgomery, and Emanuel (2012) computed drag coefficients obtained from six separate aircraft missions into the major hurricanes of Fabian and Isabel (2003), at a range of surface wind speeds from 52 to 72 m s −1 , but the data showed a large scatter (the values of C d ranged from 0.4 × 10 −3 to 5.2 × 10 −3 for wind speeds greater than 52 m s −1 ). Thus, it is hard to decide which of the candidate schemes in Figure 1 was the most realistic according to the observations available.

Simulated results
The WRF model was used to investigate how different C d parametrizations (i.e. the roughness length z 0 ) affected the two simulated super typhoons. Typhoon Haiyan (2013), which had a maximum wind speed of 87 m s −1 , and Typhoon Maysak (2015), which had a maximum wind speed of 72 m s −1 , both formed in the western North Pacific Ocean. The simulated tracks, maximum 10-m wind speeds, Charnock1955, TY2001, TYlimit, Moon2007, Davis2008, Andreas2012, Edson2013, and Liu2012 to represent the parametrizations tested. As mentioned above, the formula for z 0 in the TYlimit scheme was the same as in the TY2001 scheme except that the TYlimit scheme had an upper limit (z 0 ≤ 2.85 × 10 −3 m). The specific values of C d by the different schemes at moderate (U 10 = 30 m s −1 ) and high (U 10 = 70 m s −1 ) wind speed are presented in Table 1.
The relationships between C d and the surface wind speed in the parametrizations tested are shown in Figure 1. The value of C d was larger and increased at a higher rate in the Edson2013 scheme than in the other wind-speedrelated or friction-velocity-related schemes. C d increased Table 1. the values of C d for the different parametrizations at moderate (U 10 = 30 m s −1 ) and high (U 10 = 70 m s −1 ) wind speed.  in Figure 2(a) and (d)), which could be ignored. Compared with the observed tracks, the simulated tracks for Typhoon Maysak (shown in Figure 2(a)) were much more accurate than those for Typhoon Haiyan (shown in Figure 2(d)). All of the simulated tracks for Typhoon Haiyan were to the right of the best track, and the bias reached almost 2° of latitude after 48 h. However, the simulated tracks for Typhoon Maysak were quite close to the observed track. and minimum SLPs for the typhoons are presented here. Further structural analyses and comparisons are presented only for Typhoon Haiyan because of its extraordinary intensity.

Tracks and intensities
There were slight differences (less than 50 km) between the different simulated tracks for the typhoons (as shown the simulated typhoons reached a quasi-steady state after 54 h. RU10 and R35 were both lowest in the simulation using the Liu2012 scheme and highest in the simulation using the Davis2008 scheme.

Structure: wind and pressure
The relationships between the minimum surface central pressures and the maximum surface wind speeds for the simulations of Typhoon Haiyan and the results for the JTWC best track are shown in Figure 3(c). The simulations and observations intersected at roughly 56 m s −1 at 940 hPa, above which the wind speeds and pressures were higher in the simulations than the observations, and below which the opposite was true. The simulated points showed considerable scatter at relatively low pressures. The results in the simulation using the Moon2007, TYlimit, Andreas2012, and Liu2012 scheme were much closer to the best track, but only the simulation using the Liu2012 scheme could reproduce the points at wind speeds higher than 80 m s −1 . Four radius-height plots of the horizontal, radial, and vertical winds at 67 h are shown in Figure 3(d) as examples, allowing further assessment of the structure of the simulated Typhoon Haiyan. The maximum horizontal wind speeds in the lowest 8 km were more than 95 m s −1 in the simulations using the Edson2013, Moon2007, and Liu2012 schemes, but only reached roughly 85 m s −1 in the simulation using the Davis2008 scheme. The radii of the maximum horizontal winds were less than 50 km (approximately 40 km) to the right of the center in the simulations using the Edson2013, Moon2007, and Liu2012 schemes, and the Liu2012 scheme also gave a strong wind zone on the left side of the center. However, the radii of the maximum horizontal winds were roughly 60 km on both sides of the center in the simulation using the Davis 2008 scheme. Strong inflow and outflow occurred mainly below and above 2 km, respectively. Inflow was more obvious to the right of the center, while outflow was dominant on the left. Convergence at low levels caused upward motion near the area at which the maximum horizontal wind speed occurred, and the vertical motion was much weaker in the As shown in Figure 2(b)-(f ), both the maximum U 10 and the minimum SLP were very sensitive to the parametrization selected. Small values of C d in the simulations seemed to lead to strong surface winds. For example, the simulation using the Liu2012 scheme gave the highest maximum surface wind speed for both typhoons because it gave the lowest C d at high wind speeds. However, though the values of C d were smaller in the TYlimit and Moon2007 schemes than that in the Andreas2012 scheme at high wind speeds, the maximum surface winds from 48 to 72 h in the simulations of typhoon Maysak were weaker using the TYlimit and Moon2007 schemes than that using the Andreas2012 scheme. There were also some exceptions related to the simulated minimum SLP. For example, the Edson2013 scheme gave a larger C d but a lower minimum SLP than most of the other wind-speed-related parametrizations for both typhoons. Compared with the observations, for typhoon Haiyan, the simulation using the Liu2012 scheme gave the lowest RMSEs of the maximum U 10 and the minimum SLP (as shown in Table 2), whereas the Davis2008 scheme led the worst simulation. However, for Typhoon Maysak, the RMSE of the maximum U 10 was almost 3-5 m s −1 larger in the simulation using the Liu2012 scheme than using the other schemes, and the Liu2012 scheme also led to a bigger deviation for the SLP than most of the other schemes (except the Moon2007 scheme).

Size
The radius of the maximum U 10 (RU10) and the average radius of 35 kts (approximately 17.5 m s −1 ) winds (R35) were defined as the inner-core and the outer size of Typhoon Haiyan, respectively. The inner-core sizes of the simulations (shown in Figure 3(a)) increased slightly from 35 to 45 km on average between 12 and 72 h. However, in the simulation using the Liu2012 scheme, RU10 decreased slowly up to 50 h, and then increased slowly with time. The outer sizes (shown in Figure 3(b)) increased quickly and linearly from 150 to 325 km (350 km in the simulation using the Davis2008 scheme) between 12 and 54 h, then grew slowly to 350 km by 72 h, which suggests that   simulation using the Moon2007 scheme than in the other simulations.

Summary
Eight parametrizations of C d were comprehensively compared and assessed in this study. It is difficult to observe the behaviour of C d in strong winds, so some parametrizations of C d (e.g. Charnock 1955;Davis et al. 2008;Edson et al. 2013) have been extrapolated from field measurements in weak winds, which all reveal that C d rapidly increases as the surface wind speed increases. By combining the wave height and wave steepness, the formula proposed by Taylor and Yelland (2001) reduces the rate at which the C d increases with wind speed. Laboratory observational datasets (Donelan et al. 2004;Troitskaya et al. 2012) acquired in strong winds tend to cause C d to be saturated at wind speeds exceeding about 30 m s −1 . This feature is approximately described by the formulae proposed by Andreas, Mahrt, and Vickers (2012) and Moon et al. (2007) and the restrained formula proposed by Taylor and Yelland (2001) in the WRF model. A wave-age-related scheme (Liu, Guan, and Xie 2012) was also used in this study, and this gave a decreasing C d value at wind speeds of more than 25-33 m s −1 , consistent with field measurements (Jarosz et al. 2007;Powell, Vickery, and Reinhold 2003). The observations of C d were limited in the region and range of wind speeds, and presented big scatters, and thus more accurate observations at high wind speeds are necessary for validation of the C d parametrizations.
Numerical simulations of super typhoons Haiyan (2013) and Maysak (2015) were performed using version 3.6 of the WRF model with different C d parametrizations. Large differences between the eight schemes had almost no influence on the tracks of the simulated typhoons. The typhoon intensity indicated by the maximum 10-m wind speed and the minimum SLP was very sensitive to the selected parametrization. The size, wind-pressure relationship, and vertical structure of the wind field for the simulated Typhoon Haiyan were also affected by the surface momentum flux.
The simulations performed in this study serve to raise awareness of how and to what extent the different C d parametrizations can affect simulated typhoons, but do not reveal which parametrization can help produce the most accurate simulation compared with observation. Some other physical or dynamical parametrized processes, the spatial resolution of the WRF model, and the accuracy of initial conditions may also have a significant influence on typhoon modelling. As the next step of this research, more typhoon cases will be simulated using a coupled oceanwave-atmosphere model to give some general results.