1.  

Introduction

Heavy rainfall, strong winds and typhoon waves (hereinafter TWs) caused by a typhoon often pose a significant threat to people’s life and property near shore, as well as marine transportation and marine engineering. When a typhoon passes through the ocean, the TWs generated are highly energetic, causing devastating damages to offshore platforms, ships and coastal engineering. Studies have demonstrated an increasing trend of the potential destructive power of a typhoon under the condition of global warming (Emanuel, 2005). Therefore, timely and accurate prediction of TWs is of profound significance for disaster prevention and mitigation. Some wave models, such as WAVEWATCH III (hereinafter WW3), SWAN (the Simulating Waves Nearshore) and STWAVE (Steady-State Spectral Wave), have been widely used for the simulation and hindcast of ocean waves solely or coupled with atmospheric or oceanic component (Kim and Lee, 2018). As for a wave model, forcing wind field plays a key role, in particular, continuously transmits momentum to the ocean, driving the generation and development of TWs (Li et al., 2018; Markina et al., 2018). Cavaleri (2009) pointed that errors which were inherent in surface winds remain the main source of the uncertainty in the results of wave modelling. Therefore, an accurate representation of forcing wind fields, especially in strong wind conditions, facilitates the evaluation of the momentum flux at the air-sea interface, and thus improves the simulation and hindcast of TWs.

At present, synthetic forcing wind fields, which superimpose the wind field calculated from the parametric typhoon model on the reanalysis data, have been commonly used for the simulation of TWs and storm surges (Xu et al., 2017; Das, 2018). The parametric typhoon model, moreover, can be divided into symmetric model and asymmetric model. The idealised wind vortex taken from a symmetric wind model in particular, which was proposed by Holland (1980) (henceforth denoted H80), has been widely used as the wind forcing field. Specifically, poor approximation of the radial profile in the inner core of some storms in H80 model has been noted by Willoughby and Rahn (2004). Owing to the drift motion, the wind field is actually not symmetrical and is intensified on the right-hand side of typhoons as it should be in the northern hemisphere. Thus, some asymmetric wind field models were developed through superimposing the moving component of a typhoon to the H80 model (Jakobsen and Madsen, 2004). Gao et al., (2013) developed a revised H80 model by using a linear weighting of parameter sets computed from all available isotachs in each quadrant (i.e. generalised asymmetric Holland model, GAHM). To enable the modeling of asymmetries, the parametric radial profile of modelling inner and outer cores separately has also been explored and has shown promise to enhance overall accuracy (Jelesnianski et al., 1992; Willoughby et al., 2006; Holland et al., 2010). In addition, high-resolution H*wind (Hurricane wind) data, which was developed by NOAA-HRD (National Oceanic and Atmospheric Administration’s Hurricane Research Division, Powell et al., 1998), has also been blended with reanalysis data for building a forcing wind input for wave modelling. Its grid spacing is about 6 km and it covers the square area about 500 km away from the typhoon center. Even so, it is still not widely used as a forcing wind field because of the fact that it cannot cover the entire computational domain. Moreover, forcing wind fields derived from mesoscale atmospheric models, such as GFLD (Geophysical Fluid Dynamics Laboratory) Hurricane Prediction System (Xiang et al., 2015), MM5 (the nonhydrostatic fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model, Chen and Dudhia, 2001) and WRF (Weather Research and Forecasting) model (Bricheno et al., 2013), were also used for wave modeling solely or by coupled modeling systems (Wu et al., 2017; Zhao et al., 2017). Studies have shown that the wind output from atmospheric model can reproduce 10 m wind fields similar to observations, and the simulated wind fields have been applied for forcing spectral wave models (Langodan et al., 2014; Markina et al., 2018).

Momentum transfer at the air-sea interface also plays a very important role in the development of typhoon waves (Pianezze et al., 2018). The momentum flux in models is generally calculated by solving the wind stress coefficient (also known as momentum exchange coefficient) Cd. Earlier studies (Large and Pond, 1981; Wu, 1982) suggested a strong linear dependence of the Cd on the surface wind speed. However, in situ (Powell et al., 2003; Jarosz et al., 2007; Holthuijsen et al., 2012) and experimental measurements (Donelan et al., 2004) indicated that the empirical linear relation between Cd and surface wind speed does not apply to full wind range, especially for the strong winds. Without a better choice, a series of quadratic functions (Hwang, 2011; Zijlema et al., 2012) were proposed to fit several sets of field measured data. In order to explore the physical controls of Cd, the momentum and energy conservation equations of the turbulent stress, wave-induced stress, and viscous stress were put forward in the wave boundary layer model to obtain the momentum flux under the complex sea state (Moon et al., 2004b, hereinafter M-WBLM), but the behaviour of the Cd obtained from M-WBLM is not consistent with the measured data at medium–high wind speed. Some studies suggested that the effects of sea spray generated by breaking waves (Zhang et al., 2016) or modulation of turbulent behaviour (Chen and Yu, 2016) should be considered in the framework of M-WBLM. Based on this, the M-WBLM was extended by Reichl et al. (2014, hereinafter R-WBLM) and Chen and Yu (2016, hereinafter C-WBLM). R-WBLM modelled the spectral saturation tail level as a function of wind speed, wave age or both, and C-WBLM added the energy dissipation due to the stratification of sea spray. Both R-WBLM and C-WBLM can capture the decreasing tendency of Cd under extreme wind conditions. Fan et al. (2009) applied the M-WBLM into the WW3 model and its results yielded an improved forecast of significant wave height but an underestimation of dominant wavelength. Chen and Yu (2017) applied C-WBLM to the STWAVE model for storm wave modelling and results showed a better agreement with observed data over deep water. However, two extended WBLM have not been widely evaluated in the TWs modelling or the coupled system, especially for the R-WBLM.

Studies of TWs modelling have focused on evaluating the roles of different wind fields data (Langodan et al., 2014) or the source terms such as wind-wave interaction term, nonlinear wave-wave interaction term, and dissipation term (Seemanth et al., 2016). However, numerical simulations of the TWs are rare while both ensuring a more accurate wind input and with extended WBLM simultaneously using the WW3. Thus, we will use WW3 model to compare the performance of two constructed wind fields with one output from WRF model on the simulations of TWs. Further, we explored a modification of ST2 wind stress coefficient scheme by introducing R-WBLM to WW3. Given this, we then carried out a thorough analysis of the difference performance in the simulations of TWs between R-WBLM and the original wind stress coefficient parameterisation formula (Chalikov, 1995) within ST2. This will contribute to building a comprehensive framework for wave modelling and forecasting.

The following sections are organized as follows. Section 2 gives a brief review of Typhoon Kalmaegi and observations we collected from different platforms (buoy stations and radar altimeter). Details pertaining to the parametric wind models and WRF model configuration are also incorporated in this section, as well as the configuration of wave model and introduction of drag coefficient evaluation methods. Section 3 presents the overall analysis of the performance of wave modelling with various wind forcing fields and two different drag coefficient schemes. Section 4 provides a brief discussion and summarises the results and limitations in our study.

2.  

Data and methods

2.1.  

Kalmaegi (2014)

Severe Typhoon (STY) Kalmaegi best-track datasets are taken from Joint Typhoon Warning Center (JTWC, http://www.usno.navy.mil/JTWC) including the information of centre position, the maximum 1 min sustained wind speed and the minimum sea level pressure for each 6 h interval. Kalmaegi was generated on the Northwest Pacific Ocean and first identified as a tropical depression at 1800 UTC 10 September 2014. The maximum wind speed rose rapidly to 33 m s−1 and continued to intensify after landing the northern Philippines Island. After passing Hainan Island, Kalmaegi was classified as a Severe Typhoon with the maximum 10 min averaged wind speed of 42 m s−1 and the minimum central pressure of 960 hPa at 1800 UTC 15 September. Finally, Kalmaegi landed in Guangning province of Vietnam and died out rapidly after turning west-west north under the influence of a subtropical high-pressure cell. The best tracks of Kalmaegi and the evolution of its intensity are shown in Fig. 1. The geographical constraints include water depth and shoreline information, for which we applied the 1-min gridded bathymetry for the world (ETOPO1) (https://www.ngdc.noaa.gov/mgg/global/global.html).

Fig. 1.  

Double-layer calculation domains in WRF model (WRF-D01 and WRF-D02), single-layer calculation domain in WW3 model (WW3-D01) and terrain from ETOPO1 (unit: m); Best track and the intensity of severe typhoon Kalmaegi from JTWC, the positions of two buoys (B1 and B4) and scanning paths of SARAL (Janson-2) radar altimeter are marked by coloured dots, red triangles and black (blue) dashed lines, respectively. The intensity is classified into six categories based on maximum 10-m wind speed, i.e., tropical depression (TD), tropical storm (TS), severe tropical storm (STS), typhoon (TY), severe typhoon (STY), and super typhoon (SP-TY).

2.1.1.  

In-situ buoy measurements

There were two buoys deployed in the South China Sea (i.e. B1, located at 116 E, 19.7 N, and B4, located at 117.5 E, 19.2 N in Fig. 1). Data from buoys will be used for validating the wind forcing and the wave simulations after quality control procedures. Both of them were on the right side of the track, and the depth of the positions were 1630 m and 3690 m, respectively. The wind vector was measured every 3 min and averaged every hour with the highest value (gust) retained by a wind monitor (model 05103, R. M. Young Company, USA, http://www.youngusa.com/) from 4 m above the mean sea level, which was transformed to the 10 m height using a bulk formula (Large and Pond, 1981). Standard accuracy by this wind monitor is ± 0.3 m s−1 (0.6 mph) and ± 2 degrees for wind speed and wind direction, respectively. Significant wave height (Hs), peak wave period (Tp) and mean wave direction (Dm) from the two buoys were measured by a TRIAXYS wave sensor (AXYS Technologies, Sidney, Canada). Similarly, what we need to emphasise is that the accuracy of measuring Hs, Tp and Dm by this wave sensor can reach at 99%, 99% and 3 degrees, respectively (MacIsaac and Naeth, 2013).

2.1.2.  

Radar altimeter measurements (RAs)

It is particularly worth mentioning here that satellite altimeter missions can also provide reliable wave height measurements (Liu et al., 2016a), which have also been extensively used for validation of wave modeling (Ardhuin et al., 2010; Liu et al., 2017). In our study, we collected three geophysical records from two satellite altimeters: Jason-2 (Dumont et al., 2011) (one blue dashed line in Fig. 1) and SARAL (Bronner et al., 2013) (two black dashed lines in Fig. 1). For clarity, we rename SARAL two geophysical records to Saral-707 and Saral-722 in the light of the pass number of altimetre. Data from those two satellites are denoted as the geophysical data records (GDRs). A brief summary of these data has also been presented in Liu et al. (2016a). Table 1 summarises the detail information regarding these data, including their sources, data duration, and the orbit parameters, etc. However, before the data can be analysed, they should be carefully quality controlled, calibrated and validated. We adopt empirical criteria from Liu et al (2017, more details see their Section 2.1.2) to remove the erroneous and suspicious records. We, following Liu et al. (2016a, more details see their Table 3), use two liner corrections Eqs. (1) and (2) to calibrate Jason-2’s wave height records and SARAL’s wave height records, respectively,

Jason-2:

((1) )
Hs=1.019×Hs0.050

SARAL:

((2) )
Hs=0.997×Hs0.056

Where Hs and Hs are the raw and corrected significant wave heights, respectively.

2.2.  

Wind fields

High-quality forcing wind field is crucial for the prediction of the TWs properties. For our study, wind fields for STY Kalmaegi are developed from three different sources including symmetric wind field model, asymmetric wind field model and WRF.

2.2.1.  

Symmetric wind field model (SYM)

Holland (1980) proposed an analytical model of the radial profiles of sea level pressure and wind velocity components due to the balance of the pressure gradient, the Coriolis force, the centrifugal force and the friction force in a typhoon. The H80 model function expression is as follows:

((3) )
Vsym=Bρa(Rmaxr)B(PnPc)exp[(Rmaxr)B]+(rf2)2rf2
where Pn and Pc are the ambient undisturbed synoptic background surface pressure and the central surface pressure of the typhoon, respectively; r is the distance from the typhoon centre; ρa is the density of the air (1.15 kg m−3); f is the Coriolis force. In addition, the shape parameter B and the maximum wind speed radius Rmax are two essential parameters in the model, which control the eye diameter of the typhoon and affect the accuracy of the reconstructed wind field. Studies have applied different formulas to solve parameters B(Madsen and Jakobsen, 2004; Rego and Li, 2010). Mattocks and Forbes (2008) proposed a more detailed equation to calculate the parameter B:
((4) )
B=[(VmVT)/WPBL]2ρaePnPc
where B is limited to values between 1 and 2.5. Vm is the maximum wind speed; WPBL is a wind reduction factor whose value set to 0.9 (Powell et al., 2003); VT is the typhoon forward motion. The informations (Vm, Pc) were extracted from the JTWC best-track data. Specifically, B values are variable for entire track of when constructing wind field. The calculation results show that B values are greater than 1 but far less than 2.5 during the development and maturity of the Kalmaegi (0000 UTC 14 Sep to 1200 UTC 16 Sep). In the formation and decay stage of the typhoon, by contrast, B values are less than but very close to 1. Rmax refers to the distance between the center of a typhoon and its band of strongest winds whose value directly affects the distribution of wind speed or air pressure. Here we use the empirical equation proposed by Knaff et al. (2007, hereafter K07) based on measured data:
((5) )
Rmax=m0+m1Vm+m2(φ25°)
where m0, m1, and m2 in the equation depend on the sea area where the typhoon occurs; φ is the latitude of typhoon centre. STY Kalmaegi studied in this paper was generated and appeared in the Western Pacific Ocean. According to the rules proposed by Knaff et al. (2007), m0, m1, and m2 were set to 38.0, –0.1167 and –0.004, respectively. Pan et al. (2016) made a comparison between K07 and Xie et al. (2006, hereafter X06) for the calculation of Rmax using several typhoons. They concluded that the K07 method should be used when Vm is large, while the X06 method should be used when the Vm is small. The pity is that they did not give the specific critical value of Vm. Knaff et al. (2018) updated these three coefficients for the Western North Pacific basin. They pointed out that the original coefficients in K07 was too small and too symmetric, resulting in unrealistically small wind radii. As many researches mentioned, solving the data-driven model parameter Rmax is still an open problem (Chavas et al., 2015). Certainly, more detailed work on this issue still needs to be further studied. In order to embody H80 wind field more intuitively, we decompose wind speed Vsym into u and v components (Fig. 2). Red and blue areas indicate positive and negative extremes in the velocity field. The left (right) plot is the horizontal (vertical) mode.

Fig. 2.  

Models of H80 of Kalmaegi (1200 UTC 15 Sep). Red and blue areas indicate positive and negative extremes in the velocity field. The left (right) plot is the horizontal (vertical) mode. The black arrow shows the forward speed direction.

2.2.2.  

Asymmetric wind field model (AYM)

Due to the capability of capturing the main features of typhoon and flexibility of running at any desired resolution, the Holland model is kindly appreciated in TWs simulation. Despite the popularity of the H80, the limitation of axially symmetric about the typhoon eye is still block its practical application. To overcome this limitation, the radius of maximum wind speed (Rmax) has been replaced by an azimuthally varying Rmax(θ), where θ is the azimuthal angel around the typhoon centre (Xie et al., 2006; Mattocks and Forbes 2008). Considering the translational velocity of the typhoon centre, superimposing a storm-relative moving speed that varies with the distance from the typhoon center on the symmetric wind field model is also popularised. Lin and Chavas (2012) found that the surface manifestation of the translation vector appeared to be about 55% of the translation magnitude with direction rotated 20 degrees cyclonically in the Atlantic basin. Knaff and Zehr (2007) and Courtney and Knaff (2009) employed a slightly non-linear relationship from Schwerdt et al. (1979) that found a comparable percentage of the translation magnitude. The details of this effect can have a substantial impact on the simulation of extreme sea conditions e.g. storm surge (Lin and Chavas 2012). To remove the influence of the typhoon motion, here we use the exponential function form proposed by Jakobsen and Madsen (2004) to calculate the storm-relative moving speed:

((6) )
Vmov=Vmcexp(r/RG)
where RG is the length scale of the moving speed and is about 500 km; Vmc is the moving velocity vector of the typhoon centre whose values were calculated from best-track data about the average 6 h longitude and latitude. Here, we should point out that the relationship between Vmc and VT is VT=1.5Vmc0.63 (Mattocks and Forbes, 2008). The final expression of the asymmetric wind field can be written as:
((7) )
Vaym=Vsym+Vmov

Neither of the two analytical typhoon models can capture and reproduce the characteristics of the wind field far from the typhoon centre, but they have an effect on the improvement of accuracy near the centre of the typhoon. However, the global wind datasets, which based on satellite remote sensing assimilation including CCMP (Cross-Calibrated Multi-Platform), ECMWF (European Centre for Medium-Range Weather Forecasts) and NCEP (National Centres for Environmental Prediction), perform better in the areas farther away from the typhoon centre but poorly in reproducing the wind speed around the typhoon centre (Pan et al., 2016). So it is necessary to synthesise forcing wind input for wave simulation by blending the reanalysis data with hurricane wind model together. Compared with aforementioned reanalysis datasets, a variational analysis method was applied to the CCMP winds to assimilate extensive cross-calibrated multiple satellite datasets (e.g. SSM/I, TMI, AMSR-E, ADEOS-2 and QuikSCAT) with in situ data (e.g. ship and buoy observations) and ECMWF analyses. Even so, we still cannot neglect the shortcoming of underestimating wind events at high wind speeds (>25 m s−1) for CCMP winds (Atlas et al., 2011). Thus far, CCMP data is widely used as a forcing field for waves (Jin et al., 2014; Li et al., 2018), as well as to storm surge (Afshar-Kaveh et al., 2017) and ocean current (Li et al., 2014) regarding more available observations are integrated into it. Therefore, we use CCMP data as background winds for blended wind field as well.

Previous studies have shown that the weighting coefficient in the blended wind field is an important parameter, and its expected value is often adjusted to the size related to the maximum wind speed radius by empirical formula (Xu et al., 2007). Pan et al. (2016) determined the two critical values R1 and R2 finding a position where the difference between the maximum wind speed radius of the two wind fields was smallest, and then the superimposed wind field was obtained. Moreover, this method is flexible and also ensure a smooth transition between two wind fields. The specific expression for blended wind field Vsup can be calculated as follows:

((8) )
R1=Ropβ0Rt
((9) )
R2=Rop+(1β0)Rt
((10) )
Vsup={Vmod                 ;  r<R1(1α)Vmod+αVccmp ;  R1<r<R2Vccmp                 ;  r>R2
where β0 is the empirical parameter and was set to 0.3; Rt is the width of the transition wind field, and generally was set to 1.2Rmax; Vmod refers to Vsym or Vaym; α is a transitional function, which can be expressed as α=(rR1)/(R2R1). For detailed algorithms about Rop can be found in Pan et al. (2016). After the superposition, two modified wind fields of the typhoon Kalmaegi can be obtained. Samples of the blended symmetric wind field (blended asymmetric wind field) after superposing CCMP winds are plotted in the left column (middle column) of Fig. 4.

2.2.3.  

WRF wind field

Using the version 3.9 WRF model, a nested configuration with two domains (Fig. 1) was employed. The horizontal resolutions were 27 km and 9 km, respectively, using 44 vertical layers. The time steps used in Domain 1 (D01) and Domain 2 (D02) were 60 s and 20 s, respectively. The pressure at the top layer was 50 hPa. The initial conditions and lateral boundary conditions in the WRF model were obtained from the NCEP FNL (Final) data-set, with a 1° × 1° grid and at intervals of 6-h.

The simulation of Kalmaegi was started at 0000 UTC 12 September 2014 and ended at 0000 UTC 19 September 2014. This period covered the whole process of generation, reinforcement, landing and extinction of Kalmaegi. The hourly results from D02 are used in the analysis below. Beginning in Version 3.9, an option to use physics suites was introduced. There are currently 2 available suites (“Conus” and “Tropical”) consist of a combination of physics options. “Tropical” suite selection is based on evaluation of real-time tropical cyclone prediction over Western Pacific (2015, 2016), Eastern Pacific (2015) and Atlantic (selected storms) basins using Model for Prediction Across Scales (MPAS, Bullock et al., 2018), while “Conus” suite selection is based on many years of real-time forecasting experiments over the Continental US. Therefore, the simulation was conducted using the “Tropical” suite rather than “Conus” suite for the purpose of typhoon simulations in this study. Along with the WSM6 (WRF single-moment 6-class, Hong and Lim, 2006) and RRTMG (Rapid Radiative Transfer Model for GCMs, Bae et al., 2016) schemes for the microphysics and radiation simulations, the physical options New Tiedtke (Zhang and Wang, 2018) cumulus parameterisation scheme, the Yonsei University Planetary Boundary Layer (PBL) with nonlocal turbulent mixing (Hong et al., 2006), the old MM5 Monin–Obukhov surface layer scheme (Salesky and Chamecki, 2012) and the Noah land surface model (Chen and Dudhia, 2001) in the “Tropical” suite were used. In addition, sea surface temperature from Real-Time Global (RTG) with 0.5° horizontal grid resolution was integrated. We will focus our analysis on the results of the simulation for D02.

The simulated tracks of Kalmaegi driven by the “Tropical” suite selection is shown in Fig. 3a along with the observed best track from JTWC data. At the beginning of the initial stage of the simulation (first 48 h), the simulated tracks for STY Kalmaegi is in good agreement with the observed track, and the bias reached minimal (roughly 23.81 km) at the time instant of 48 h. After 96 h integration corresponding the first landing time instant (0000 UTC 16 September), the simulated positions of typhoon centre gradually lagged behind the observed positions, but still moved along the observed trajectory. Since then, the bias had been larger reaching almost 2.89°. The results of minimum sea level pressure (Pmin) and maximum 10-m wind speed (U10, max) are shown in Fig. 3b and c, respectively. The simulated Pmin is higher than the observed value at the initial period. With the increase of typhoon intensity, the results of simulated Pmin and observed data intersected at running roughly 72 h. By contrast, the results of simulated Pmin are higher than the observations before 72 h, and it gradually begin to be lower than the observed value and the typhoon intensity are overestimated after 72 h. Figure 3c shows that most of the time the results of simulated U10, max are higher than the observed value, which demonstrates a good correspondence with the trend of Pmin from Fig. 3b. According to statistical analysis on the data, the average RMSE of U10, max and Pmin are 7.36 m s−1 and 4.12 hPa, respectively.

Fig. 3.  

(a) Observed best track of Kalmaegi from JTWC (red line) and the simulated track using ‘Tropical’ suite (blue line). (b) The minimum sea surface pressure, with lines as defined in (a). (c) As in (b) but for the maximum 10 m wind speed.

Since the numerical results of the “Tropical” suite, to some extent, are in good agreement with the observations with relatively small intensity and track errors, the wind field output from the “Tropical” suite selection will be used as the third forcing wind for TWs simulation.

2.2.4.  

Comparisons of three wind fields

In order to visualise the spatial difference of three wind fields more intuitively, three time instants (0000 UTC 15 September, 0600 UTC 15 September and 1200 UTC 15 September) are selected for qualitative comparisons of three wind fields about spatial structure (shown in Fig. 4). The wind fields constructed by the symmetric model (Fig. 4a–c) depict an obvious symmetrical spatial structure at different moments. The asymmetry at three time instants (Fig. 4d–f) are expected to be induced by typhoon’s forward motion, resulting in a high wind speed zone at the front right side of the system in the forward direction. In contrast, the maximum wind speed of WRF output (Fig. 4g–i) from these areas can be captured on the right side or right rear side of typhoon centre. Meanwhile, the differences of asymmetrical spatial structure between the asymmetric model wind field and WRF output are obvious, in particular, the localised distribution of the maximum wind speed zone at 0000 UTC 15 September and 0600 UTC 15 September. The spatial structure of the three wind fields are quite distinctive, and the numerical simulation results for the TWs modelling driven by different wind fields are bound to be different.

Fig. 4.  

Comparison of spatial distributions among symmetric wind field (left column: a–c), asymmetric wind field (middle column: d–f), and WRF output (right column: g–i) at 0000 UTC 15 Sep (top row: a, d, g), 0600 UTC 15 Sep (middle row: b, e, h) and 1200 UTC 15 Sep (bottom row: c, f, i). Buoys location are marked by black triangular and the center of the typhoon are marked by black circle.

Wind data at B1 was invalid due to some unknown reasons, so the validation was tested only by B4. Wind speed produced from three wind field and observed wind speed, together with CCMP at station B4 are shown in Fig. 5. As can be seen from Fig. 5, when the typhoon centre is far from the B4, all of the wind speed data are quite consistent with the observed wind speed. However, when the typhoon centre is very close to the B4, the differences of the peak wind speed among different wind fields, which the maximum wind speed of observation is up to around 27 m s−1, are particularly obvious. To be sure, CCMP data does underestimate peak wind speed while the other three different wind fields improve the wind speed at different levels during this period. Wind speeds obtained from asymmetric wind field are closer to the values of observed compared with that from symmetric wind speed owing to typhoon’s forward motion in the asymmetric wind model. Ideally, the peak wind speed from WRF’s output is rather close to the measurement at B4 since the “Tropical” suite was used for the WRF modeling and a quiet satisfactory result was presented with the correlation being up to 0.96 (confidence interval is 5%).

Fig. 5.  

Time series of 10 m wind speed from CCMP and three wind fields versus observed data at B4.

2.3.  

Wave modelling

2.3.1.  

Numerical wave model

In our study, WW3 (V5.16) is used for wave modelling. WW3 is a full spectral third-generation wave model that solves the balance equation for the wave action density spectrum in the conservative. Several physical processes are generally considered in net source terms (S) including atmosphere-wave interaction (Sin), a nonlinear wave-wave interaction (Snl) and a wave-ocean interaction term that generally contains the dissipation (Sds). Sin and Sds represent separate processes, but are often interrelated owing to the balance of these two source terms governs the integral growth characteristics of the wave energy (Tolman, 2016). Several source term packages (Sin+Sds) are available including WAM cycle 3 (Komen et al., 1984), TC96 (Tolman and Chalikov, 1996), WAM4+ (Bidlot et al., 2005), A10 (Ardhuin et al., 2010), and BYDRZ (e.g. Babanin et al., 2007, 2010) which can be employed in model by activation of ST1, ST2, ST3, ST4 and ST6 switches, respectively. For the detailed description of each package, the reader is referred to υσερ μανυαλ ανδ ρεϕερενχεσ τηερειν. In order to facilitate the following research (see Section 2.3.3), we will select ST2 only to launch discussing.

2.3.2.  

Experimental designs

We used a single domain in the wave model, between longitudes 90°–165° E and latitudes 0°–45° N as illustrated the red rectangle frame in Fig. 1. The spatial domain was discretised using a rectangular grid with resolution of 0.125° × 0.125°. The frequency range was 0.0285 to 1.1726 Hz with an increment factor of 1.1 and 72 directions with 5° resolution. The maximum global time step, the maximum CFL (Courant – Friedrichs – Lewy) time step for xy and k-theta and the minimum source time step were defined as 1200 s, 400 s, 1200 s and 300 s, respectively. The simulation is initialized at 0000 UTC 1 September 2014 and 11 days ahead for model spinup. In order to compare the performance of numerical wave modeling when using different forcing wind fields, four experiments (test0, test1, test2 and test3, see Tab. 2) are carried out firstly.

2.3.3.  

Drag coefficient evaluation methods

In general, the calculation of wind stress vector τ is considered as to determine the wind stress coefficient Cd or friction velocity u due to these parameters are related to each other through:

((11) )
τ=ρaCd|U10|U10;  τ=|τ|=ρau2
where ρa is the density of the air; U10 is the wind speed vector at a 10 m height above the sea surface; τ is the size of τ. In WW3 model, different source term packages use different formulas to calculate u or Cd. ST1 uses Wu’s formula (1982); ST2 uses algorithm proposed by Chalikov (1995); the Quasi-liner theory is used both in ST3 and ST4 packages (Janssen, 1991); Hwang’s (2011) scheme is used in ST6 package. In our study, in order to investigate the application of R-WBLM under strong wind condition, one of source term packages need to be chosen to carry out a contrastive study between R-WBLM and self-contained formula in WW3.

In particular, we would like to highlight that the R-WBLM subroutines have been added into the WW3 (V5.16). The provided subroutines allow the computation of the total wind stress based on the WW3 wavenumber-direction spectrum. Nevertheless, R-WBLM subroutines are provided for use as part of wave-ocean or other coupled systems. It’s designed to compute additional quantities (e.g. surface roughness length, Charnock parameter, friction velocity) related to the surface wave field which are intended to be passed to external models (e.g. ocean models). The motivation for these subroutines is to allow the external model to include the impact of waves on quantities (see Section 2.5.1 of V5.16 manual). However, we run WW3 separately without the second model at present, so R-WBLM needs to be embedded into one of the source term packages based on its dependence on the wave spectrum firstly.

The input peak frequency fp,i defined as Eq. (12) is a vital parameter in R-WBLM, which is not available without a wave spectrum. It will directly affect the calculation of cutoff wave number (Moon et al., 2004b) and a further impact on the various stresses within the R-WBLM. In ST2 source term package, this frequency is estimated from the equivalent frequency of the positive part of the input source term as Eq. (12). It is difficult to get an estimate value from the remaining four packages (ST1/3/4/6) due to the complex relevance among parameters within the wind input term.

((12) )
fp, i=f2cg1max[0,Swind(k,θ)]dfdθf3cg1max[0,Swind(k,θ)]dfdθ
where Swind is the wind energy input source term for wave energy conservation equation, cg is group velocity. For above purpose, R-WBLM eventually was embedded into the ST2 package. The computation formulas of wind stress coefficient about R-WBLM and Chalikov’s (1995) are described as following section.

2.3.3.1.  

R-WBLM

The transfer of momentum is generally equivalent to the wind stress τ and expressed as

((13) )
τ=τt(z)+τw(z)
where τt and τw are the local turbulent stress and wave-induced stress, respectively. The wind profile within the wave boundary layer is calculated from the kinetic energy conservation equation:
((14) )
ddz(Uτ)+dΠwdz+dΠtdzρaε=0
where U is the mean wind vector; Πw and Πt are the vertical transport of energy related to wave-induced motion and turbulent motion, respectively; ε is the viscous dissipation rate. The main idea of R-WBLM to solve the wind stress is to enter the 10 m wind velocity above the sea surface U10, in for the first step, and then use the wave spectrum generated from WW3 and an initial estimation of viscous stress τυ to calculate wind profile and predict U10, out by integrating. The iteration is required to continue until the final value of the wind speed converges, that is U10,inU10, out, and the wind stress τ within the wave boundary layer is finally determined. Finally, according to the Eq. (11), the wind stress coefficient is obtained and then provided to WW3 for further wave field calculation.

2.3.3.2.  

Chalikov’s formula

((15) )
Cd=103[0.021+10.4/(R1.23+1.85)]

where

((16) )
R=max{0,ln[z10g/(χαU102)]}
where χ=0.2 is a constant, α=0.57(max(u,0.0001)/cp)1.5 (Janssen, 1989) is the conventional nondimensional energy level at high frequencies.u and cp are friction velocity and spectrum peak phase speed. According to Eq. (11), Cd is finally obtained by several times of integrating.

For simplicity, the R-WBLM embedded into ST2 is recorded as ST2+. ST2+ and ST2 will be used to represent the same source term packages (Sin+Sds) but with different wind stress calculation methods. It should be noted that the level of the saturation spectrum in R-WBLM has a significant effect on the values of drag coefficient. In our study, we set this tuning parameter via namelist FLD1 as TAILLEV = 0.006. The remaining three tuning parameters, which are also crucial to the results of drag coefficient, name ‘High Frequency Tail Method’, ‘Tail transition ratio 1’ and ‘Tail transition ratio 2’, respectively. The values of them were set to 1, 1.250, and 3.000, respectively. To further investigate the performance of employing ST2+ into WW3, test4 are added in Tab. 2 (see Section 3.2).

3.  

Results

3.1.  

Comparisons of forcing wind fields

Firstly, in order to obtain a better quantitative and qualitative perspective of three wind fields performance for wave parameters (Hs, Tp,Dm) while using ST2, detailed analyses were performed throughout comparisons between measurements and simulated results as illustrated in Fig. 6. The whole period of typhoon process was divides into three time segments. One thing that needs to be explained here is that the mentioned time segments defined as the time when STY Kalmaegi was close to, near, and away from two buoys. The first time segment is from at 0000 UTC on 12 Sep to 1300 UTC 14 Sep (short for Before-TY); the second time segment is from 1300 UTC on 14 Sep to 0000 UTC on 16 Sep (short for During-TY); the third time segment is from 0000 UTC on 16 Sep to 0000 UTC on 19 Sep (short for After-TY). Moreover, it has to be kept in mind that the contribution of CCMP data in each time segment was markedly different, so we added a controlled experiment (test0 in Tab. 2) as mentioned in Section 2.3.2. Three statistical parameters: Standard Deviation std, RMSEε0, and correlation coefficient ρ, are utilised to evaluate the performance of the different forcing wind fields:

((17) )
std=1N1iN(xix¯)2
((18) )
ε0=1Ni=1N(xiyi)2
((19) )
ρ=i=1N(xix¯)(yiy¯)i=1N(xix¯)2i=1N(yiy¯)2
where x and y represent the simulated and measured wave quantities, the bar over x and y denotes their mean value, and N is the sample size of the collocations. Statistical results for B4 and B1 stations are presented in Tabs. 3 and 4, respectively.

Fig. 6.  

Time series of simulated Hs (m) (left column), Tp (s) (middle column), and Dm (°) (right column) vs. the observed data at: B4 (top row) and B1 (bottom row).

During the period of “During-TY”, statistical results show that the simulated three wave parameters driven by WRF’s output with measurements at B4 have the smallest ε0 (0.913 m) and highest ρ (0.875). It is also clearly that the application of WRF’s output results in the smallest ε0 (0.817 m, 11.857°) and highest ρ (0.858, 0.983) for Tp and Dm at B4. Similarity, the smallest ε0 (1.169 m) and highest ρ (0.945) for Hs are obtained from WRF’s output at B1. During the period of “Before-TY” and “After-TY”, it can be seen that application of all wind fields lead to a general agreement between simulation and observation. In terms of the numerical simulation results of the Hs during “After-TY” period, to our surprise, the results show that: symmetric wind field is the best, the second is asymmetric wind field, the worst is WRF’s output. We know that different physics schemes (e.g. microphysics, longwave/shortwave radiation and cumulus parameterisation) in WRF model can lead to different results. The poor performance using WRF’s output during “After-TY” period is a sign that WRF modelling need to be further employed to choose a better set of physical options for a better-quality wind forcing, but we do not want to put too much concentrate on the sensitivity of intensity and track to these physics schemes.

Now we are paying our attention to test0’s results. We note that a noticeable underestimation of Hs and Tp driven by CCMP data during the period of “During-TY”, as seen from left and middle column in Fig. 6, which is mainly attributed to the deficiency of tending to underestimate high winds. Unlike time segment “During-TY”, the differences with two other time segments in simulating wave parameters are not very obvious (see Tabs. 3 and 4). Combined with Fig. 5, we can conclude that CCMP can provide high-quality wind data while wind speed are less around 15 m s−1. This feature is clearly reflected from Fig. 6. Specifically, it is obvious that the experimental results generally are in good agreement with buoy observations during periods of “Before-TY” and “After-TY”.

Eventually, the preliminary conclusion is that the best performance wind field is WRF’s output during “During-TY” among these cases. In the following sections, we will use WRF wind forcing for our next research based on measurements from various platforms.

3.2.  

Drag coefficient Cd

The variations of the wind stress coefficient Cd against the wind speed resulted from ST2+ and ST2 at all space points within the simulated time period and several data set of field observations (Powell et al., 2003; Donelan et al., 2004; Jarosz et al., 2007; Holthuijsen et al., 2012) are shown in Fig. 7. It is clear that ST2 scheme significantly overestimates Cd when compared with measurements and with an upper limit of 0.0058 for wind speed more than 27 m s−1. In contrast, the variation of Cd calculated by ST2+ is very close to the observation results and show a good decreasing tendency of Cd at medium–high wind speed.

Fig. 7.  

Wind stress coefficient obtained from ST2+ case (red dots), ST2 case (black dots) and field observations (e.g. Powell (with error bars); Jarosz; Donelan; Holthuijsen (with error bars)) vary with 10 m wind speed.

3.3.  

Calibration against measurements

According to the data acquisition moment at Kalmaegi’s lifetime from RAs (see the fourth column records in Tab. 1), model results at 2200 UTC 14 Sep, 1000 UTC 15 Sep and 2100 UTC 15 Sep these three moments have been used for the convenience of study. The comparison of spatial distribution of the simulated Hs obtained from ST2+ and ST2 schemes at these three representative time instants as shown in Fig. 8. Subgraphs (g–i) are the differences of Hs between ST2+ and ST2 scheme at the same time instant, with 10 m wind filed contours and scanning paths of RAs are superimposed. All maximal values of Hs appear in the right forward quadrant along the typhoon track (as indicated in Fig. 8a–f), and this phenomenon is again demonstrated that higher waves are more likely to appear on the right hand of the typhoon track, because the waves in this region are subjected to longer forcing of strong wind. Again, attention is drawn to an obvious difference of the spatial distribution of Hs between ST2+ and ST2. We can conduct a further preliminary conclusion that the wind stress coefficient calculation schemes will have an effect on the numerical simulation of wave heights, especially for the high wind sea areas.

Fig. 8.  

Spatial distribution of Hs driven by WRF’s output among ST2+ case (top row: a–c), ST2 case (middle row: d–f) and differences between two cases (bottom row: g–i) at 2200 UTC 14 Sep (left column a, d, g), 1000 UTC 15 Sep (middle column b, e, h) and 2100 UTC 15 Sep (right column c, f, i). The heavy solid lines denote the track of typhoon. 10 m wind speed contours are showed with grey lines in g, h, i (contour intervals of the wind field are 10 m s−1). Bold dashed lines from left column to right column denote scanning paths from Saral-707 (black), Saral-722 (black), and Jason-2 (blue) at corresponding moment, respectively.

Scanning path from Jason-2 just passed through the centre of the typhoon (see Fig. 8c, 8f), while the other two scanning paths from SARAL appeared before (see Fig. 8a, 8d) or after (see Fig. 8b, 8e) the centre of the typhoon. It is obvious that ST2+ scheme (top row) yields smaller Hs than that from ST2 scheme (middle row) at the area between 20 and 30 m s−1 wind speed contours, which are located on the right hand side of the typhoon track. This is evidently because ST2+ yields smaller Cd at these areas when the wind speed is around 20–40 m s−1 (see Fig. 7). The results show that the differences of simulated Hs between ST2+ and ST2 (|Hs,ST2+Hs,ST2|) are up to 1 m, or even higher at 1000 UTC 15 Sep and 2100 UTC 15 Sep.

Comparisons about Hs between RAs and model results from ST2+ and ST2 are illustrated in Fig. 9a–c. Sixty-two and 98 collocations are obtained from Saral-707 and Saral-722, while 91 collocations are given by Janson-2. Combined with RMSE ε0 and correlation coefficient ρ, another two statistical parameters bias b and scatter indexSI are also utilised to quantify model skill:

((20) )
b=1Ni=1N(xiyi)
((21) )
SI=1Ni=1N[(xix¯)(yiy¯)]2y¯

Fig. 9.  

Comparison between the calibration results for the Hs* measurements according to Eqs. (1–2) and the simulated results from ST2+ (blue dots) and ST2 (red dots), respectively. (a–c) Saral-707, Saral-722 and Jason2, respectively. Similarity, comparison between buoy-based measurements and the simulated results are showed in (d) and (e). Error statistics for bias b,RMSE ε0, correlation coefficient ρ, and scatter index SI are given inset. The dashed line is 1:1 line.

Error statistics are given inset in Fig. 9. In general, agreements between numerical results from two schemes and all three altimeters are good with high correlation (0.74) and a small scatter index (0.27). It is easy to detect that simulated Hs obtained from ST2+ show the lowest RMSE (0.66 m, 1.02 m and 1.69 m) and scatter index (0.11, 0.19 and 0.26) compared with that from ST2 (Fig. 9a–c). For Fig. 9a, simulated Hs obtained from ST2+ displays a negative bias value (–0.55) compared with Saral-707-calibrated Hs from collocations while results obtained from ST2 shows a smaller negative bias value (–0.39). Although, both of them show the same values in correlation coefficient and scatter index. This indicates that ST2+ and ST2 provide similar performance in computing wave heights when Hs < 6 m. Most importantly, the performance of ST2+ scheme is relatively encouraging under the high sea states (Hs > 6 m). Specifically, inspection of Fig. 9c shows that the agreement between simulated Hs obtained from ST2+ and Jason2-calibrated Hs measurements is better than the results from ST2 scheme when Hs > 6 m. This conclusion, to some extent, do also hold for Fig. 9b. The results here indicate that the ST2+ yields smaller Cd than ST2 and therefore should be expected to yield reduced wave heights under high wind speed conditions.

Unlike RAs, buoy-based measurements obtained from a duration of time on fixed positions, so high sea states were less frequently collected by those two buoys (see Fig. 9d–e). Even so, agreement between model results from ST2+ and buoys is slightly better than the results from ST2. Among the four statistics considered in Fig. 9d–e, model results from ST2+ have slightly smaller biases (0.01 m and 0.05 for B4 and B1, respectively) than that from ST2; while other three error statistics show a similar feature. To summarise briefly, based on the comparisons between model results and measurements from various platforms, we concluded that ST2+ show more accuracy than ST2 in specification of wave height parameters.

3.4.  

Spatial distribution of wave age

The differences of simulated Hs between ST2+ and ST2 at high wind speed areas is may attributed to the effect of the sea spray generated by breaking waves under extreme sea conditions. Much of what distinguishes ST2+ from ST2 is that it included the impact of sea state by considering different growth rates for waves. In related studies, the effect of swell was simply ignored and the growth rate coefficient cβ was set as 40 (Moon et al., 2004b) or 25 (Chen and Yu, 2016). However, cβ is modified to explicitly account for the swell in R-WBLM as follows:

((22) )
cβ={25                         :cos(θθw) ;c/u*l<1010+15cos[π(c/u*10)/15]:;10<c/u*l<255:;25<c/u*l25:cos(θθw)<0
where ul=(τt/ρa)0.5 is the local turbulent friction velocity; θ and θw are wind stress direction and wave direction, respectively. It was classified into three regimes, wind forced waves (c/ul<10), the swell forcing wind (c/ul25) and a transition regime (10c/ul<25). Thus, while the wind stress based parameterisation of ST2 prevents swell from having an outstanding impact on the wind stress, ST2+ introduces a significant contribution from swell to the wind stress.

Observational evidences show that the wind stress may also be modulated by the presence of swell (García-Nava et al., 2012). Figure 10 gives the spatial distribution of wave age simulated by ST2+ (a) and ST2 (b) at the same time step with that in Fig. 8c. The wave age β=cp/U10 is used to describe the growth state of the wave, and cp is the phase velocity of the waves at the spectral peak frequency. It is generally defined that the wave age exceeds 1.2 for swells, while less than 1.2 for wind forced waves. The white region in Fig. 10 is equivalent to β=1.2. Both two buoys were deployed in central South China Sea and waves generated on the Pacific Ocean was blocked by the northern Luzon and the group of small islands located between Taiwan and Luzon. So the effects of swells were relatively limited for two buoys. In addition, wave states near buoys are wind forced waves with young wave age in the process of typhoon center approaching buoys, because the waves in this region are subjected to longer forcing of strong wind. This conclusion, to some extent, do also hold for Jason2 measurements (see Fig. 8c). It is clearly that the sea states are dominated by wind forced waves within the scope of area shown in this figure. From this point of view, R-WBLM results show more sea state dependence while the sea state dependence of ST2 is much weaker. According to the previous studies, sea spray generated by breaking waves which almost always generated locally lead to energy dissipation and make relatively less momentum exchange at air-sea interface (Moon et al., 2004a; Chen and Yu, 2017), so we can have an intuitive view that the areas controlled by wind-wave component have more effected by R-WBLM. To summarise, there is no doubt that the R-WBLM takes better account of the sea state dependence in terms of calculating wind stress.

Fig. 10.  

Spatial distribution of wave age (coloured) between ST2+ (a) and ST2 (b) cases at 2100 UTC 15 Sep, respectively. The heavy solid lines denote the track of typhoon. Wind speed contours are showed with grey lines, black ‘+’ and black ‘●’ denote the position of U10,max and Hs,max at this moment, respectively.

4.  

Discussion and conclusions

In this study, the application ability of different wind field models and wave boundary layer model (WBLM) are investigated simultaneously in WW3 model using typhoon Kalmaegi. As controlling vital factor for typhoon waves hindcast, a better-quality forcing wind field was chosen from constructed typhoon models and mesoscale atmospheric model (WRF), and also a modification of drag coefficient scheme in ST2 source term package was done by introducing R-WBLM to WW3 for TWs simulation. Severe typhoon Kalmaegi (2014) was used as a test case since observations from various platforms, including two buoys measurements and three radar altimetre records, were available for validation. The main findings can be drawn as follows:

  1. Two constructed wind fields based on parametric model are constructed and one from the WRF’s output using “Tropical” suite were validated. In addition to providing insight into two reconstructed wind forcing fields based on H80, the present analysis also provides a comprehensive description of the wind output from mesoscale atmospheric models WRF. Especially, physics suite “Tropical” has been provided in version 3.9 since we always confuse the combination of physics options when configuring the WRF model for a typhoon hindcast. Such a combination scheme of parameters has important consequences for WRF modelling. The results from WRF modelling show that the simulated track and intensity of typhoon present a reasonable numerical results. We use measurement of wind speed at B4 to validate WRF modelling results and a quiet satisfactory result is obtained with the correlation being up to 0.96 (confidence interval is 5%). Hence, the wind field from WRF output is also reasonable as a forcing wind input. As a result, three wind fields improve the local wind speed in various degrees but with different spatial distribution, especially for the location of maximum wind speed.
  2. It was not surprising that the underestimation of Hs and Tp driven by CCMP data during the period of “During-TY” probably resulted from the deficiency of tending to underestimate high winds for CCMP data. However, CCMP can provide high-quality wind data while wind speed are less around 15 m s−1. Such results are generally consistent with the reported by Jin et al. (2014) and Li et al. (2018). We can see clearly the model results are in good agreement with buoy observations in our experiments (see Fig. 5). Based on the comparisons between model results (e.g. Hs) driven by three wind fields and two buoy-based measurements, we concluded that simulated Hs driven by WRF’ output are more consistent with measurements during the time of “During-TY”: for B4 and B1, the correlations are 0.875 and 0.945, respectively. Hence, wind from WRF’s output was used for the further experiments, which mainly focused on the application of different wind stress coefficient schemes.
  3. Based on the requirement an estimation for the peak frequency fp, i in internal framework of R-WBLM, R-WBLM subroutines were embedded into ST2. Figure 7 gives the variations of the wind stress coefficient Cd against the wind speed resulted from R-WBLM and Chalikov’s formula. It is obvious that the variance of Cd in the range of medium and high wind speed is unreasonable. Therefore, ST2 package is not recommended when using WW3 to simulate TWs in the future research. The decrement of saturation level with increasing wind speed, which was considered into the framework of WBLM by Reichl et al. (2014), contribute to reducing the Cd at high wind speed. Hence, there is no doubt that the R-WBLM is recommend to calculate wind stress coefficient, as the Cd obtained from it is qualitatively consistent with many sets of field observations under extreme sea conditions.
  4. Comparisons about Hs between measurements from calibrated RAs and buoy-based and model results from ST2+ and ST2 were implemented. We collected wave height measurements from radar altimeters: Jason-2 and SARAL to enrich validation data. It has also been demonstrated that Radar altimetre measurements can also provide high-quality wave height measurements by Zieger et al (2009) and Liu et al (2016a). The performance of ST2+ scheme is relatively encouraging under the high sea states ( Hs > 6 m) (this can be found from error statistics given by Fig. 9). Specifically, agreements between simulated Hs obtained from ST2+ and Saral-722, Jason2-calibrated Hs measurements are better than the results from ST2 scheme when Hs>6 m (see Fig. 9b,c). The results reconfirm that the ST2+ yields smaller Cd than ST2 and therefore should be expected to yield reduced wave heights under high wind speed conditions.

Even so, there are still many limitations in aspect of verifying the superiority of applying R-WBLM in our numerical simulation case. Naturally, the accuracy of wind forcing still needs to be improved, which is associated with physical options in WRF modelling (Li et al., 2018; Zhang and Wang, 2018). In addition, the superiority of R-WBLM in this case cannot be verified more quantitatively due to limited observations data. It is still difficult to conclude that the WBLM is a better choice for wave modelling based on the numerical results for only a case occurred in the Northwest Pacific basin. To generalise the findings presented above, more well-observed sample of typhoons should be simulated and analysed. This is especially necessary for the operational wave forecasting purpose (Chao and Tolman, 2010; Liu et al., 2017). Although the application effect of WBLM in numerical simulation of TWs in the Northwest Pacific basin is hard to be verified by more real cases, such a dominant role of WBLM for TWs modelling has previously been reported by Chen and Yu (2017) and Du et al. (2019). Similarly, they applied WBLM to storm wave modelling by taking cases of Hurricane Katrina (2005) occurred the Gulf of Mexico and two North Sea storms, respectively. Results from their research concluded that the WBLM is conclusively shown to be a valuable option for storm wave modeling.

Having stated the findings of the present study, it is undeniable that the WBLM has increasingly become a recognised method for evaluation of the wind stress under extreme sea conditions. Generally speaking, R-WBLM not only embodies a rather good performance in the aspect of variations of the wind stress coefficient against the wind speed, but also takes better account of the sea state dependence for TWs simulation at the air-sea interface. In order to demonstrate the improvement of the numerical simulation when using R-WBLM, more cases need to be verified in the future. As an important factor in the air-sea interaction process, the momentum exchange is of great significance to the evaluation and improvement of the simulation precision of typhoon track and intensity, even for upper ocean mixing and cooling, as well as the improvement of the numerical simulation of TWs. Therefore, it is necessary to further investigate and verify the application of R-WBLM in wave model or even in atmosphere-wave-ocean coupled models.