Effects of Drag Coefficients on Surface Heat Flux during Typhoon Kalmaegi (2014)

Considering the increasing damages made by the landfalls of tropical cyclones (TCs) and the large uncertainties of TC intensity in the warming ocean, it is important to accurately forecast the TC track and intensity. However, compared with the shrinking errors in TC track forecasts, predicted TC intensity lacks striking improvements over recent decades, which could be partly attributed to the limited understanding of energy exchange in the form of surface heat flux (i.e., sensible and latent heat fluxes)(Emanuel, 1986; Jaimes et al., 2015) between TCs and the ocean. Given that the surface heat flux is the primary energy source for TC intensification (Barnes and Dolling, 2013), it is closely correlated with the surface drag coefficient ($ {C}_{\mathrm{d}} $), which determines not only the momentum transport but also thermal and moisture fluxes through the exchange coefficients ($ {C}_{\mathrm{h}} $ and $ {C}_{\mathrm{q}} $) over the air–sea interface. Unfortunately, lacking knowledge of the processes between TCs and the ocean poses an obstacle to fully understanding the impact of $ {C}_{\mathrm{d}} $ on surface heat flux. Within this context, how exactly $ {C}_{\mathrm{d}} $ behaves under strong cyclonic winds and what the effects are on surface heat flux and thus the development of TCs remain open questions.

Due to the destructive winds and dense clouds of TCs, effective observations of complex air–sea interactions in the boundary layer are scarce. Therefore, $ {C}_{\mathrm{d}} $ in high wind speed scenarios (>22 m s⁻¹) is usually measured by extrapolation from low or medium wind speeds (Black et al., 2007). With GPS dropwindsonde data, Powell et al. (2003) revealed that $ {C}_{\mathrm{d}} $ became uniform and then slightly decreased at wind speeds greater than 33 m s⁻¹. A similar conclusion was drawn by Donelan et al. (2004) with laboratory extreme wind experiments. Shay and Jacob (2006) further confirmed a “saturation” $ {C}_{\mathrm{d}} $ at 30 m s⁻¹ by equating internal wave–ocean fluxes to surface winds at the 10-m level. In addition, Vickery et al. (2009) utilized GPS dropwindsonde data and found that $ {C}_{\mathrm{d}} $ appeared to level off and then decrease at wind speeds greater than 60 m s⁻¹. Most recently, Zou et al. (2018) used two subsurface buoys during Typhoon Megi (2010) to estimate $ {C}_{\mathrm{d}} $ at high wind speeds through two bottom-up methods based on the turbulence closure scheme and bulk model. They found that $ {C}_{\mathrm{d}} $ increased to a maximum at a critical wind speed of approximately 30 m s⁻¹ and then became uniform with wind speed. Although there are no common views of the threshold value at which $ {C}_{\mathrm{d}} $ begins to decrease, it is widely accepted that $ {C}_{\mathrm{d}} $ could not increase monotonically with increasing wind speed, which has been parameterized in the numerical models (Bell et al., 2012).

Presently, the uncertainties of $ {C}_{\mathrm{d}} $ at high wind speeds support the use of three different possible parameterization schemes in numerical models. The first and most widely used $ {C}_{\mathrm{d}} $ trend is monotonically increasing with wind speed based on the Charnock relationship (Charnock, 1955), though some studies have demonstrated that this trend is not suitable in strong TC scenarios (Green and Zhang, 2013; Zou et al., 2018). The second scheme, which is more consistent with observations (Powell et al., 2003; Donelan et al., 2004), maintains a constant sea surface roughness under high winds up to 40 m s⁻¹. Such a scheme is currently predominant in TC simulations, where wind speeds reach up to 70 m s⁻¹ (Drennan et al., 2007; Jeong et al., 2012; Richter and Stern, 2014). The third scheme produces a reduced $ {C}_{\mathrm{d}} $ at high wind speeds. Considering the effects of the wavy ocean surface, Liu et al. (2011) implemented the sea-spray-affected aerodynamic roughness ($ {z}_{0} $, which is equivalent to $ {C}_{\mathrm{d}} $) in a coupled atmosphere–wave–ocean model and further confirmed that TCs weakened with the increasing $ {C}_{\mathrm{d}} $ at low-to-medium wind speeds but intensified in high-wind conditions when $ {C}_{\mathrm{d}} $ decreased. Taking the feedbacks of sea surface temperature (SST) and near-surface wind into consideration, whether this scheme adequately represents $ {C}_{\mathrm{d}} $ in TC conditions remains to be investigated.

Limited by insufficient in situ observations and the uncertainties in $ {C}_{\mathrm{d}} $ under high-wind conditions, investigations on the sensitivity of surface heat flux or TC intensity to $ {C}_{\mathrm{d}} $ are still in their nascent stages (Montgomery et al., 2010; Green and Zhang, 2013; Smith et al., 2014). Using an axisymmetric numerical model with balanced boundary-layer dynamics, Emanuel (1986) presented a highly influential perspective wherein the TC intensity is proportional to the square root of the ratio of surface exchange coefficients of enthalpy, $ {C}_{\mathrm{k}} $, and momentum, $ {C}_{\mathrm{d}} $. With the first direct measurements of enthalpy flux in a hurricane boundary layer, Zhang et al. (2008) found that there was no statistically significant dependence of these bulk exchange coefficients on wind speed. They concluded that the enthalpy flux into the hurricane boundary layer, which is required to initiate and sustain hurricane development, may have to come from sources other than air–sea turbulent fluxes. Based on their conclusion, the Emanuel model assumptions should be revisited. Smith and Montgomery (2010) suggested that the increase of $ {C}_{\mathrm{d}} $ around an inertial regime leads to enhanced radial inflow and, thus, significantly intensified surface heat transfer. Zhang et al. (2015) also suggested that enhancing surface friction would enhance boundary layer inflow, which in turn strengthens the convergence of angular momentum and intensifies a TC. However, this could not lead to a positive relationship between $ {C}_{\mathrm{d}} $and the intensification rate of a TC since surface friction also enhances the momentum and heat dissipation to boundary layer winds, thus dynamically causing a negative effect on TC intensification (Stern et al., 2015; Heng and Wang, 2016; Wang and Heng, 2016). Montgomery et al. (2010) suggested a $ {C}_{\mathrm{d}} $ of 2.0 $ \times $ 10⁻³ to be the threshold value determining the TC intensification rate and mature intensity in a cloud-resolved numerical model. With a $ {C}_{\mathrm{d}} $ of less than 2.0 $ \times $10⁻³, both intensification rate and mature intensity of the simulated TC increased slightly with increasing $ {C}_{\mathrm{d}} $, but the mature intensity decreased with either larger or zero $ {C}_{\mathrm{d}} $. However, some previous studies argued that both simulated TC intensification rate and mature intensity were insensitive to randomly perturbed $ {C}_{\mathrm{d}} $ with changes of over 60% (Thomsen et al., 2014; Peng et al., 2018). Chen et al. (2018) suggested the effect of $ {C}_{\mathrm{d}} $ on TC evolution was negligible because $ {C}_{\mathrm{d}} $ induced moderate sea-surface cooling (SSC) despite the dominant role of $ {C}_{\mathrm{h}} $ and $ {C}_{\mathrm{q}} $ in surface heat fluxes. Therefore, the effects of the uncertainties in $ {C}_{\mathrm{d}} $on both surface heat fluxes and TC intensity need to be better understood.

Observations from five cross-shaped arrays of buoys and moorings in the northern South China Sea (SCS, Fig. 1) precisely captured the intensity change concerning the ocean responses and associated air–sea interactions of typhoon Kalmaegi (2014) (Zhang et al., 2018). Kalmaegi crossed the northeastern end of Luzon and entered the SCS at 1800 UTC 14 September with a speed of 30 km h⁻¹, heading towards Hainan Island. It passed right over the array of buoys and moorings from 0000 UTC 15 September to 0000 UTC 16 September and thus provided us the opportunity to study the distribution and variation of heat flux under typhoon conditions. In tandem with the high-resolution air–sea coupled model, a series of sensitivity experiments are conducted to examine the effects of $ {C}_{\mathrm{d}} $ on surface heat flux distribution and variation and the feedbacks on typhoon evolution. The rest of the paper is organized as follows: The data, detailed model configurations, and experimental designs used in this study are introduced in section 2. The model output is verified with both conventional and buoy observations in section 3. The sensitivity of the surface heat flux to both SSC and $ {C}_{\mathrm{d}} $ is analyzed in section 4. Finally, conclusions and a brief discussion are given in section 5.

Figure 1.  Comparisons between the best track data and simulated results of Kalmaegi (2014). The positions of the typhoon at 0000 UTC on each date are marked by circles and squares in (a), while five buoys/moorings are numbered and marked by triangles. (b) shows the time series of the maximum 10-m wind speed of Kalmaegi (2014) from 0600 UTC 12 to 0000 UTC 17 September.

Two sets of 6-hourly best-track data for Kalmaegi are obtained from the Joint Typhoon Warning Center (JTWC) and Shanghai Typhoon Institute, China Meteorological Administration (STI/CMA), including the location (in longitude and latitude) of the TC center and the maximum sustained 10-m wind speed. The daily SST data at 9-km resolution are obtained from the optimum interpolated SST (OI_SST) of Remote Sensing Systems.

The in situ observations are from an array containing five cross-shaped buoys and four subsurface moorings (Fig. 1). Stations 1 and 4 are located to the right of Typhoon Kalmaegi (2014), while stations 2 and 5 are to the left. A detailed description of the observational design can be found in Zhang et al. (2016). Since the temperature data at station 1 was lost and the buoy at station 3 was adrift during the passage of Kalmaegi around 0000 UTC 15 September, analyses are mainly based on the data of stations 2, 4, and 5. The in situ observations include: (1) 10-m wind speed, 2-m air temperature, and SST to verify the coupled model outputs; (2) subsurface temperature down to 40 m to estimate the ocean heat content and the thermal structure in the upper layer of the ocean; and (3) near-surface humidity to calculate the latent heat flux based on three $ {C}_{\mathrm{d}} $ schemes.

The Weather Research and Forecasting (version 3.7.1 of WRF_ARW, hereinafter referred to as WRF) model (Skamarock et al., 2008) and the Stony Brook parallel ocean model (sbPOM) (Jordi and Wang, 2012), representing the atmospheric and oceanic components, respectively, were coupled to each other through the Model Coupling Toolkit (MCT) (Jacob et al., 2005; Larson et al., 2005). As a fully compressible, non-hydrostatic primitive equation atmospheric model in terrain-following vertical coordinates, WRF is widely used in regional mesoscale atmospheric modeling and operational numerical weather forecasts. The physical parameterization schemes in WRF are the same as those in (Zhang et al., 2019). WRF has 51 vertical half-sigma levels to the height of 50 hPa. There are 570 $ \times $ 652 grid points in the x and y directions, respectively, with grid spacing of 9 km, and the outermost domain covered the northwest Pacific (1°−49°N, 100°–149°E, in Fig. A1.). The atmospheric initial and lateral boundary fields adopt the 6-hourly National Centers for Environmental Prediction (NCEP) Final (FNL) operational global analysis data on 1-degree by 1-degree grids.

Figure A1.  Configuration of the coupled model domain. The color shades denote land and ocean topography in meters.

In the coupled model, ocean feedbacks are from the Stony Brook Parallel Ocean Model (sbPOM), a paralleled primitive equations ocean model with free surface and sigma coordinate. Vertically, 40 full-sigma levels are configured with a maximum depth of 5000 m. Therein, 16 levels are set in the upper 200 m to better resolve the ocean mixed layer. To reduce the computational cost and avoid unnecessary interpolation error, sbPOM is set with the same domain and horizontal resolution as the WRF model. The initial and lateral boundary conditions of temperature, salinity, currents, and sea level are acquired from the Hybrid Coordinate Ocean Model (HYCOM) + The Navy Coupled Ocean Data Assimilation (NCODA) system Global 1^o/12 analysis, which assimilates the satellite altimeter observations and in situ SST, as well as the in situ vertical temperature and salinity profiles from the Expendable bathythermographs (XBTs), ARGO floats, and moored buoys.

During the coupled simulations, the time steps of atmospheric and oceanic components are set to 30 s and 180 s, respectively. The WRF and sbPOM models interact every 1800 s, which means that the WRF model is integrated over every 60 steps to make it coincide with the sbPOM model after every 10 integration steps. The atmosphere component transfers atmospheric forces (including longwave and shortwave radiation flux, sensible heat flux, latent heat flux, rainfall, and wind stress) to the ocean. Simultaneously, SST and current feedbacks are transferred to the atmospheric component to adjust the bottom conditions and the wind-current shear (Liu et al., 2013). The specific model settings are listed in Table 2.

Theoretically, $ {C}_{\mathrm{d}} $ is a function of wind speed, sea surface state, and atmospheric stability, defined as:

where $ k $ is the von Kármán constant and $ {z}_{\mathrm{r}\mathrm{e}\mathrm{f}} $ is the reference height (often 10 m). It not only impacts the amount of momentum flux transferred into water columns but also determines the heat-flux exchange coefficients in air–sea interactions (Donelan et al., 2004):

where $ {z}_{\mathrm{h}} $ and $ {z}_{\mathrm{q}} $ are the thermal roughness length and the moisture roughness length, respectively.

To involve all the $ {C}_{\mathrm{d}} $ behaviors mentioned above, sensitivity experiments are designed using three main schemes (Fig. 2). The first two $ {C}_{\mathrm{d}} $ parameterization schemes are chosen to represent the increasing and leveling off $ {C}_{\mathrm{d}} $ with TC wind speed, which are listed in the WRF name-list file as the parameter isftcflx= 0 and 1, respectively. Option 0 is inherited from the Mesoscale Model 5 (MM5) model and yields a monotonic increase of $ {C}_{\mathrm{d}} $ with wind speed (Charnock, 1955). It is expressed as a function of $ {z}_{0} $ and has been widely used in earlier simulations of TCs with the MM5 model (hereinafter referred to as Charnock55):

Figure 2.  Plots of (a) aerodynamic roughness $ {z}_{0} $, (b) exchange coefficient of drag $ {C}_{\mathrm{d}} $，sensible heat $ {C}_{\mathrm{h}} $ and water vapor $ {C}_{\mathrm{q}} $, and (c) exchange coefficient ratio $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $ associated with 10-m wind speed for S0, S1 and S2.

where $ {u}_{*} $ is frictional velocity and $ g $ is the acceleration due to gravity. This scheme is called S0, indicating a rougher sea surface with increasing wind speed.

When the parameter isftcflx is set as 1 (the green line in Fig. 2), $ {C}_{\mathrm{d}} $ increases with increasing wind speed but saturates at the speed of 33 m s⁻¹, which is more consistent with observations (Powell et al., 2003; Donelan et al., 2004). The momentum roughness length $ {z}_{0} $ is defined as:

where

This scheme is called S1, and all other settings are the same as S0.

Considering the effect of sea spray in high-wind conditions, the third scheme used in this study was proposed by Liu and Guan (2007) (hereinafter referred to as Liu07 and shown by the blue line in Fig. 2). In this scheme, when 10-m wind speed is below 25 m s^–1, $ {z}_{0} $ and $ {C}_{\mathrm{d}} $ increase steadily with increasing wind speed. Subsequently, $ {z}_{0} $ gradually reaches its peak value at a wind speed of around 33 m s^–1. From thereon, $ {z}_{0} $ (same tendency as $ {C}_{\mathrm{d}} $) becomes uniform or decreases with increasing wind speed. This scheme is given by:

where $ {\alpha }_{*}={c}_{\mathrm{p}}/{u}_{*} $, $  \omega =\text{min}\left(1,{\alpha }_{\mathrm{c}\mathrm{r}}/\kappa {u}_{*}\right) $, and $ {\alpha }_{\mathrm{c}\mathrm{r}}=0.64\;{\rm{m}}\,{\rm{s}}^{-1}$. $ {\alpha }_{*} $ denotes the wave age, and $ {c}_{\mathrm{p}} $ denotes the heat capacity at constant pressure for dry air. This scheme is called S2.

Denoted by a black dashed line in Fig. 2a, both thermal and moisture roughness length ($ {z}_{\mathrm{h}} $ and $ {z}_{\mathrm{q}} $, respectively) are set to 10^–4 to guarantee the independent effects of $ {C}_{\mathrm{d}} $ on changes of both sensible heat $ {C}_{\mathrm{h}} $ and latent heat $ {C}_{\mathrm{q}} $ fluxes [Eqs. (2) and (3)]; i. e., $ {C}_{\mathrm{h}} $ equates to $ {C}_{\mathrm{q}} $ in the same $ {C}_{\mathrm{d}} $ scheme. To compare the relative effects of $ {C}_{\mathrm{d}} $ and SSC on surface heat fluxes, six experiments using three $ {C}_{\mathrm{d}} $ schemes respectively are designed and classified into two groups (Table 1). One group, using the atmosphere-only model, is identified with the letter U. Another group is identified with the letter C, denoting the air–sea coupled experiments. All experiments share the same domains and the same atmospheric conditions, except for the $ {C}_{\mathrm{d}} $ scheme and ocean responses. Experiments US0 and CS0 adopt the S0 $ {C}_{\mathrm{d}} $ scheme, while US1 and CS1 are conducted with the S1 $ {C}_{\mathrm{d}} $ scheme, correspondingly. These two groups are designed in order to compare their results with those of past studies (Chen et al., 2018). To analyze the effects of the spray-affected sea surface on the heat fluxes (Fig. 2), experiments using spray-involving scheme S2 are conducted with both WRF-only and coupled models and are named US2 and CS2, respectively. All the simulations are integrated for 120 hours from 0000 UTC 12 September to 0000 UTC 17 September, covering the period during which Typhoon Kalmaegi (2014) quickly developed and passed through five cross-shaped buoy and mooring arrays (Fig. 1a).

Kalmaegi (2014) formed in favorable environmental conditions over the Northwest Pacific Ocean and developed into a tropical storm on 12 September. It then took a northwestward trajectory to Luzon Island and became a Category 1 typhoon on 14 September, peaking at an intensity of 36 m s⁻¹ (Fig. 1a). Though Kalmaegi (2014) weaken0.8ed slightly during its first landfall over the Philippines, it quickly regained an intensity of 38 m s⁻¹ after entering the SCS on 15 September and began heading toward the southeast coast of China. As a fast-moving typhoon (defined as having a translation speed of > 5.6 m s⁻¹ at a latitude between 15°N and 20°N ) (Holland, 1993; Lin et al., 2009), Kalmaegi (2014) passed the mooring arrays with a translation speed of above 8.3 m s⁻¹, with its maximum wind velocity increasing gradually to around 42 m s⁻¹ in the early morning of 16 September. Regarding the maximum 10-m wind speed (VMAX, Fig. 1b), some discrepancies exist during this period between the two best track data sets. The VMAX values from CMA maintained 42 m s⁻¹ after a quick intensification and were generally stronger than the VMAX values from JTWC. After making landfall in Guangxi Province, Kalmaegi (2014) decayed quickly (not shown).

The recorded tracks of Typhoon Kalmaegi (2014) from JTWC and CMA are largely overlapping, especially when Kalmaegi (2014) enters the SCS (Fig. 1a). Like the best track data, all the simulated tracks pass the array of mooring buoys with a relatively slower translation speed than that seen in the observations, but this does not affect our findings. The overlapping tracks also suggest that the track of Kalmaegi (2014) is mainly determined by the steering flow and environmental systems and is not sensitive to SSC and $ {C}_{\mathrm{d}} $ (Green and Zhang, 2013; Chen et al., 2018).

Similar to what is seen in the best track data, the VMAX results of the six experiments reveal two different stages during the lifetime of Kalmaegi (2014) (Fig. 1b). Before 1800 UTC 14 September, all simulated TCs intensify steadily until Typhoon Kalmaegi (2014) makes its first landfall. After entering the SCS (around 0000 UTC 15 September), Kalmaegi (2014) resumes intensifying in all simulations, reaching its maximum intensity at 0000 UTC 16 September. In general, the coupled model reproduces weaker maximum intensities than the WRF-only simulations, except for those using the S1 scheme; CS1 and US1 simulate similar wind speeds (shown as the green lines in Fig.1b). Compared with other experiments, the coupled simulation with the Charnock55 scheme (CS0) produces the proximal intensity to the best tracks.

To further validate the coupled model and interpret the effect of $ {C}_{\mathrm{d}} $ on ocean response, SSTs from the coupled simulations are compared with the WindSat SSTs (Wentz et al., 2013) shown in Fig. 3. Evident cooling appears to the right of the track, with a minimum value of 26°C in the WindSat data, which is 1°C–2°C colder than the surrounding waters. Though all coupled simulations capture SSC to the right of the track, both magnitude and coverage of the SSC are weaker than the observation. This could also be a reason for a simulated TC that is stronger than the best track data suggest. Discrepancies also exist between the three coupled experiments, in which SSC in CS0 is more evident and consistent with the WindSat data, indicating stronger ocean responses with increasing $ {C}_{\mathrm{d}} $.

Figure 3.  SST (shaded) at 0000 UTC 17 September and typhoon track (black solid lines) from (a) CS0, (b) CS1, and (c) CS2, and (d) SST (shaded) at 0000 UTC 17 September from Windsat and typhoon tracks from CMA (thick solid line) and JTWC (thin solid line), respectively. Black triangles denote the positions of the five buoys/moorings.

The simulated 10-m wind speed and differences between SST and 2-m air temperature are compared with the buoy observations in Figs. 4a1–a3 and Figs. 4b1–b3, respectively. To reduce the instantaneous error in the comparisons, mean values from the nine model grid cells nearest to the buoy location are adopted. Since the sizes of the simulated TCs (represented by the radius of maximum 10-m wind speed) are similar (not shown), their effects are not considered here. For 15 September, the coupled model reproduces a tendency of the 10-m wind similar to that of the buoy observations, though in relatively higher values. To the right of the track, the buoy at station 4 captured a VMAX of about 22 m s^–1 at 1200 UTC 15 September, which is much stronger than that at stations 2 and 5. The peak value at station 2 occurred at 0000 UTC 16 September because of its further west location and the intensifying typhoon. Comparing all simulations, the 10-m wind speed of CS0 is the weakest due to the increasing $ {C}_{\mathrm{d}} $.

Figure 4.  Time series of 10-m wind speed (a1–a3), temperature differences between 2 m and sea surface (b1–b3), and latent heat flux (c1–c3) from the coupled simulation cases (CS0–2) and observational data (Obs) at stations 2, 4, and 5 from 0000 UTC 15 September to 0000 UTC 17 September. For (c1–c3), bold lines indicate coupled-model results (CS0–2), and buoy observations calculated with the three $ {C}_{\mathrm{d}} $ schemes (S0–2, named Obs_S0, Obs_S1, and Obs_S2, respectively) are shown in light-colored solid lines.

For the temperature differences ($ \Delta T $) between 2-m height and the sea surface, relatively large discrepancies exist among model results and buoy observations. At stations 2 and 5, the modeled $ \Delta T $s are similar and steady over time with a minimum value of –1°C, whereas the observed $ \Delta T $ decreased rapidly to –2.5°C on 15 September and varied significantly afterwards. This suggests that the near-surface air cooling, mostly induced by heat dissipations, evaporation, and heavy rainfall under the eyewall, could possibly be underestimated in the numerical model, especially when SSC is indistinctive under a fast-moving typhoon. However, $ \Delta T $ tendencies at station 4 were different from those at stations 2 and 5, which were to the left of the track. Positive $ \Delta T $ values of about 1.5°C were observed throughout this stage. The observed values are about 1.0°C higher than the results of experiments CS0 and CS2 and are even higher than the results of experiment CS1. This could be attributed to the underestimation of the SSC relative to the observations when the wind speed is weaker than 30 m s⁻¹ (Fig. 3), especially when using the leveling off $ {C}_{\mathrm{d}} $ scheme. In other words, SSC plays the key role in determining $ \Delta T $ under typhoon conditions on the right side of the TC but is most likely underestimated by the air–sea coupled model, even with the increasing $ {C}_{\mathrm{d}} $ scheme.

Latent heat flux (LHF), known as the primary energy source for TCs, is also investigated (Figs. 4c1–c3). LHF is written as:

where $ \rho $ is the air density and $ {L}_{\mathrm{v}} $ denotes the specific heat capacity of air at constant pressure. $\varDelta q=\left({q}_{\mathrm{s}}-{q}_{\mathrm{a}}\right)$, where $ {q}_{\mathrm{a}} $ is the specific humidity at a reference height of 2 m above the sea surface and $ {q}_{\mathrm{s}} $ is the saturation humidity at the sea surface. $ {U}_{\mathrm{a}} $ is the wind velocity at the same level of $ {q}_{\mathrm{a}} $. Since the buoys could not measure LHF directly, the changes of which was calculated with 10-m wind speed and $ \delta $q from the buoy observations according to Eq. (10) while $ {C}_{\mathrm{q}} $ estimated from the parameterization schemes. As shown in Figs. 4c1–c3, the calculated LHFs with the three $ {C}_{\mathrm{q}} $ schemes are almost overlapping, indicating the insensitivity of LHF to $ {C}_{\mathrm{q}} $. The modeled LHFs exhibit synchronous changes but are higher than the observations, especially at stations 2 and 5, which is possibly due to the relatively stronger 10-m wind speeds (Figs. 4a1–a3) and higher $ {q}_{\mathrm{a}} $ over the warmer sea surface. The modeled LHF at station 4 reaches a maximum value that is 100 W m⁻¹ higher and about 9 h earlier than the buoy observations and then quickly decreases after 1800 UTC 15 September. The LHF of CS0 shows more consistencies with the observations due to the stronger SSC compared to the other experiments. The above analyses also suggest that the $ {C}_{\mathrm{d}} $-induced ocean responses associated with LHF could be different on both sides of the typhoon track. Supported by earlier work, “right bias” is a typical phenomenon for fast-moving TCs in which wind vectors turn clockwise and strengthen the wind velocity to the right of the track. Therefore, the ocean responses could be enlarged due to the uncertainty of $ {C}_{\mathrm{d}} $, resulting in evident discrepancies in LHF.

Determining both momentum (e.g., wind stress) and enthalpy flux through the air–sea interface, $ {C}_{\mathrm{d}} $ not only directly changes near-surface wind speed, but also drives the mixing and upwelling processes in the ocean upper layer and thus the SSC under TCs. Both SSC and wind speed are known as the main reasons for changing surface heat flux. Therefore, the relative effects of the $ {C}_{\mathrm{d}} $-induced SSC and wind speed on the surface heat flux should be examined.

Since CS1 and US1 simulate similar wind speeds (Fig. 1b) for 15 September, the daily-averaged sensible and latent heat fluxes from these two experiments are firstly compared with the NCEP data to verify the effects of SSC (Fig. 5). Although the heat fluxes from the NCEP data are characterized by homogeneous distributions due to the relatively low resolution (1° × 1°), both magnitude and pattern of the data are more consistent with the results from CS1. US1, with fixed SST, overestimates the heat fluxes, especially around the five mooring buoys, further emphasizing the importance of ocean feedbacks in TC simulation (Elsberry et al., 1976; Price, 1981; Shay et al., 1992).

Figure 5.  Daily averaged sensible (upper row) and latent heat fluxes (bottom row) from (a, d) CS1, (b, e) US1, and (c, f) NCEP (shaded, units: W m⁻²) on 15 September. Black triangles denote the positions of the five buoys/moorings.

Figure 6 shows the differences in sensible and latent heat fluxes (S/LHD), SST, and 10-m wind speed between US1 and CS1. The S/LHD are characterized by negative values, especially to the right of the track, which is consistent with the remarkable SST cooling of –0.6°C to –1.0°C and slightly weakened wind speed (< –1 m s⁻¹) in CS1. Positive S/LHD are located where SSC is smaller than –0.4°C, showing an indistinctive relationship with the 10-m wind speed (about 0.4 m s⁻¹, not shown). Note that the variation of LHF is three times larger than that of SHF, although their patterns are similar. LHF is changed significantly by SSC (Wada et al., 2018), which determines the intensification of Typhoon Kalmaegi (2014) in the following stages.

Figure 6.  Daily averaged (a) sensible heat flux differences (shaded, units: W m⁻²) with SST differences (blue lines, units: °C) and (b) latent heat flux differences (shaded, units: W m⁻²) with 10-m wind speed differences (black lines, units: m s⁻¹) between US1 and CS1 on 15 September. Black triangles are the positions of the five buoys/moorings. Red squares are the positions of typhoon centers in CS1.

Figure 7 compares the different effects from $ {C}_{\mathrm{d}} $-induced SSC and wind speed on the LHD at stations 2, 4, and 5. Shown by the red lines, the differences between the SST-fixed experiments US0 and US1, US1 and US2, respectively, denote the effect of $ {C}_{\mathrm{d}} $ on LHD through wind changes. Similarly, comparisons of C0 minus C1 and C1 minus C2, respectively, reveal the effect of different $ {C}_{\mathrm{d}} $ schemes on LHD. Three green lines show the difference of LHD between the coupled simulations and the uncoupled ones, including C0 (CS0 minus US0), C1 (CS1 minus US1), and C2 (CS2 minus US2), respectively, indicating the impacts of $ {C}_{\mathrm{d}} $-induced SSC on LHD. At station 4, significant LHD (shown by green lines in Fig. 7) appears at around 1200 UTC 15 September with a magnitude much larger than what could be the effects from $ {C}_{\mathrm{d}} $-induced wind speed and different $ {C}_{\mathrm{d}} $ schemes, suggesting that the SSC, relating to the ocean status, is the predominant factor determining the air–sea surface heat fluxes. The changes of LHF through $ {C}_{\mathrm{d}} $-impacted wind speed is almost equivalent to those induced by different $ {C}_{\mathrm{d}} $ schemes, indicating that the near-surface wind speed only has a significant effect when SSC is inconspicuous. The results of Jaimes et al. (2015) also support the finding that the enthalpy fluxes are closely related to the upper oceanic thermal structure but independent of wind speed, which could be altered by the pressure field when SSC occurs. Compared with station 4, station 2 was located in the front left of the typhoon. Except for the strengthening SSC under the typhoon, cold water to the right rear of Kalmaegi (2014) was transported cyclonically to station 2, resulting in obvious LHD 12 hours later than at station 4. Located to the rear left of the typhoon, SSC at station 5 was not significant, and the LHD induced by SSC and wind speed was therefore of equal magnitude.

Figure 7.  Time series of differences in latent heat fluxes (LHD) from C_d-induced wind speed (US0–US1; US1–US2), SSC (C0=CS0–US0; C1=CS1–US1; C2=CS2–US2), and SSC change (C0–C1; C1–C2) on stations 2, 4, and 5 from 0000 UTC 12 September to 0000 UTC 17 September.

As shown by the blue lines in Fig. 7, the LHD induced by the three $ {C}_{\mathrm{d}} $ parameterization schemes in the coupled model are of similar magnitude as the wind-impacted LHD (shown by red lines), suggesting there is less of an effect by different $ {C}_{\mathrm{d}} $ schemes on enthalpy flux than by wind-driven evident SSC. To further understand the sensitivity of enthalpy flux under the storm to the three $ {C}_{\mathrm{d}} $ schemes, differences in daily surface heat flux between CS experiments are compared in Fig. 8.

Figure 8.  Daily averaged differences of (a, b) sensible and (c, d) latent heat fluxes (shaded, units: W m⁻²) and SST (blue lines) and 10-m wind speed (black lines, units: m s⁻¹) between CS0 and CS1 (left column) and CS1 and CS2 (right column) on 15 September. Black triangles indicate five mooring buoys and red squares are typhoon centers.

During the passage of typhoon Kalmaegi on 15 September, SST differs slightly between the coupled simulations (less than –0.2°C) to the right of the track, whereas the difference of 10-m wind speed maintains –1 m s⁻¹ to –3 m s⁻¹ (Figs. 8c and d), at a similar magnitude as that in Fig. 6b. Different from the pattern seen in Fig. 6, the weakened 10-m wind overlaps negative S/LHD zones randomly on both sides of the track rather than to the right of it, supporting the findings of the wind effects on LHF with weak SSC suggested in section 4.1.

Compared with Figs. 8b and d, the negative S/LHD zones of CS0 and CS1 (Figs. 8a and c) are more notable, which could be explained by the differences of $ {C}_{\mathrm{d}} $ values and their impacts during this period. Before 1200 UTC 15 September, the VMAX of CS0 (Fig. 1b) accelerates to approximately 33 m s⁻¹, with $ {C}_{\mathrm{d}} $ being around 3 × 10⁻³ (Fig. 2). Simultaneously, $ {C}_{\mathrm{d}} $ (or $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $) in CS1 gradually rises but is limited to 2.5 × 10⁻³ (or $ {C}_{\mathrm{h}}/{C}_{\mathrm{d}} $ at 0.7). The difference of $ {C}_{\mathrm{d}} $ (or $ {C}_{\mathrm{h}}/{C}_{\mathrm{d}} $) between CS0 and CS1 thus enlarges with the increasing 10-m wind speed (Figs. 2b and c) to about 0.9 × 10⁻³ at 0000 UTC on 16 September. In contrast, the difference of $ {C}_{\mathrm{d}} $ (or $ {C}_{\mathrm{h}}/{C}_{\mathrm{d}} $) between CS1 and CS2 decreases first with increasing wind speed and then increases moderately to 0.6 × 10⁻³, being 0.3 × 10⁻³ smaller than the difference of $ {C}_{\mathrm{d}} $ between CS0 and CS1. Since the increased surface friction could decelerate the near-surface wind speed immediately and thus inhibit the enthalpy flux, which is much faster than the ocean cooling, smaller S/LHD show up, suggesting high sensitivity of surface heat fluxes to $ {C}_{\mathrm{d}} $ schemes.

To further compare the sensitivity of LHF to $ {C}_{\mathrm{d}} $ and SSC in the inner core of Kalmaegi , LHD between CS0 and CS1, and CS1 and US1, respectively, are normalized by the azimuthally averaged LHF of CS1 over the period from 0000 UTC 15 September to 0000 UTC 17 September (Fig. 9). Following Jaimes et al. (2015), in this study, the inner core of a TC is defined as four times the radius of maximum winds (4Rmax).

Figure 9.  Azimuth–time Hovmöller diagrams of the normalized differences of latent heat flux between (a) CS0 and CS1 and (b) CS1 and US1 (shaded, units: W m⁻²). Normalization is done by the azimuthally averaged latent heat fluxes in CS1 from Rmax to 4Rmax over the period 0000 UTC 15 September to 0000 UTC 17 September. Note the different label bars of the two snapshots.

During the early period when the VMAX is still weaker than 30 m s⁻¹, positive LHD between CS0 and CS1 show up in front of the TC (Fig. 9a). As Kalmaegi (2014) intensifies, $ {C}_{\mathrm{d}} $ increases quickly with increasing wind speed, especially for CS0. LHD between CS0 and CS1 exhibit negative values in most areas east of the track, supporting the above findings that LHF relates closely to both wind speed and $ {C}_{\mathrm{d}} $ when SST changes moderately. In contrast, evident SSC dominates LHD between CS1 and US1 (Fig. 9b), which are negative in all directions, though both TCs are of similar intensity, especially to the northeast of the center. Moreover, LHD induced by $ {C}_{\mathrm{d}} $ parameterization scheme differences (Fig. 9a) are one-third or less than LHD caused by SSC (Fig. 9b), suggesting the wind-induced SSC is the more dominant factor for LHF compared to wind speed.

The downward momentum transport through an air–sea interface, determined by $ {C}_{\mathrm{d}} $ and wind speed, usually induces cooling of the ocean mixed layer (OML) when the SST becomes 0.5°C colder than the water below (e.g., Black et al. 2007) and vertical mixing and upwelling of deeper water occurs. Following Zhang et al. (2016), sea temperature of the upper 40 m is extracted from all coupled simulations and compared with the buoy observation at station 4. Figures 10a-d present the temperature differences between run time and 0000 UTC 15 September. Time series of temperature in the ocean upper 160 m are preseted in Fig. A2 . Before 0600 UTC 15 September, the simulated sea temperature below 20 m increases by about 0.5°C, while the surface layer cools slightly, indicating weak vertical mixing in the shallow OML. Later on, evident cooling occurs in the OML as the typhoon intensifies. As shown by Fig. 2b, although the wind speed of CS0 is weakest at station 4,$ {C}_{\mathrm{d}} $ is largest there throughout 15 September. Contrarily, $ {C}_{\mathrm{d}} $ in CS1 is the smallest at station 4, with the wind speed being weaker than 33 m s⁻¹, within the same period. The cooling OML (Figs. 10a–c) is proportionate to $ {C}_{\mathrm{d}} $ on 15 September for both CS0 and CS1, and CS0 produces the most significant cooling with the longest duration, while CS1 exhibits the weakest effects. Compared with the OML pattern at station 4, the magnitude of the cooling in CS0 is most consistent with the buoy observation. Note that the observed maximum cooling occurred at the depth of 20 m, different from the sea surface in a coupled model. This could be attributed to the incorrectly resolved solar heating in the model or horizontal diffusion caused by mesoscale vortexes and highlights the urgent need in developing an advanced air–sea coupled model. The variation of ocean heat capacity (OHC) with time could further reflect the effect of $ {C}_{\mathrm{d}} $ on the thermal status of the OML (Fig. 10e) (Jaimes et al., 2015). The model outputs a similar tendency with, but relatively higher than, the observed $ \varDelta $OHC, which could be caused by the existing cold eddies observed by the altimeter at station 4 (not shown). In CS0, OHC decreases gradually by 100 KJ cm⁻² after 0600 UTC 15 September, while the OHC of CS1 simultaneously decreases only by about 50 KJ cm⁻².

Figure 10.  Time series of temperature differences in the ocean upper 40 m from CS0 (a), CS1 (b), CS2 (c), and buoy observations (Obs) at station 4 (d) and ocean heat capacity (OHC) differences (e) between 0000 UTC 17 September and 0000 UTC 15 September.

To further elucidate the effects of $ {C}_{\mathrm{d}} $ on the slow recovery ocean responses and the surface heat fluxes evolution, the same analyses shown in Fig. 8 are made for 16 September and shown in Fig. 11. Since the maximum intensity of Kalmaegi (2014) on 0000 UTC 16 September varies greatly between the three experiments, the differences of SST, wind field, and enthalpy flux between the CSs are thus significant near 114°E, 19°N. The notable positive differences of SST, combined with the heat flux differences between CS1 and CS2 to the right of the center (Fig. 11b) are related to the quickly decreasing $ {C}_{\mathrm{d}} $ with peak wind speed of higher than 40 m s⁻¹. Across the mooring array, however, the daily averaged 10-m wind speeds differ little, while the SST differences among the three experiments remain larger than 0.3°C (Figs. 11a and b) and are well correlated with the variations of heat fluxes. Therefore, the distribution of surface heat flux could be greatly affected by choice of $ {C}_{\mathrm{d}} $ parameterization scheme through the longer lasting SSC, even after the passage of a TC.

Figure 11.  Same as Fig. 8, but for 16 September.

Choosing a suitable $ {C}_{\mathrm{d}} $ parameterization scheme is of great importance for numerical models, as $ {C}_{\mathrm{d}} $ not only controls the surface wind stress directly but also affects the thermal status of the upper ocean layers. With observed $ {C}_{\mathrm{d}} $ and $ {C}_{\mathrm{k}} $ profiles being incorporated in numerical models, better TC intensity and structure could be simulated, especially over the ocean (Ming and Zhang, 2016). However, the uncertainties of $ {C}_{\mathrm{d}} $ in high-wind scenarios make the estimation of surface heat fluxes in numerical simulations challenging, limiting improvement of typhoon intensity forecasting. Benefiting from buoy observations (Zhang et al., 2016) and the fully-coupled WRF-sbPOM model, a series of $ {C}_{\mathrm{d}} $-sensitivity numerical simulations were conducted in this study to understand the impacts of $ {C}_{\mathrm{d}} $ on surface heat fluxes during Typhoon Kalmaegi (2014).

The air–sea coupled model was first confirmed to perform better in simulating the intensity of Kalmaegi with more reasonable SSC than the SST-fixed WRF model by comparing output with two best track datasets and the satellite-observed SST. Sensitivity experiments with three $ {C}_{\mathrm{d}} $ parameterization schemes produced evident differences in near-surface winds and air–sea heat flux, suggesting the $ {C}_{\mathrm{d}} $ scheme to be an influential factor impacting numerical simulations. Due to the significant SSC associated with an evident reduction in heat flux, the coupled model with the increasing $ {C}_{\mathrm{d}} $ parameterization scheme reproduced the weakest TC, which was most consistent with the best track data. The reason is twofold: 1) as the coupled model tended to overestimate the intensity of a Category 1 typhoon, the increasing $ {C}_{\mathrm{d}} $ scheme could effectively enhance the surface friction and SSC, therefore reducing the intensity; 2) the ratio of $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $ in the CS0 experiment was between 0.5 and 0.65, which had been calculated with the first direct measurements of enthalpy flux in the hurricane boundary layer (Zhang et al., 2008) and proven to be suitable in estimating the air–sea fluxes relative to the OHC changes (Jaimes et al., 2015). Nevertheless, the $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $ ratios in CS1 and CS2 were greater than 0.7 as wind accelerated, especially between 15 September and 16 September, when the $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $ ratio exceeded 0.8 in CS2. It was significantly greater than the steady value of 0.65 in Zhang et al. (2008). With $ {C}_{\mathrm{d}} $ leveling off or decreasing in the high-wind scenario, SSC was inhibited. OML cooling and heat loss through the air–sea interface was therefore weakened, resulting in a stronger TC.

Compared with the buoy observations in the SCS, the simulated 10-m wind speeds exhibited a similar tendency but at a relatively higher magnitude than the buoy observations. The $ \Delta T $ between 2-m air temperature and SST was also compared between observations and model outputs. Discrepancies existed on different sides of the TC track. To the left of Kalmaegi (2014), observed $ \Delta $T declined more obviously than the model results, suggesting that the near-surface air cooling could be likely underestimated in the numerical models. Rather, the significantly higher $ \Delta T $ observed at station 4 indicated the underestimation of SSC to the right of TC track by the coupled model. Due to the warmer initial SST in the WRF-sbPOM model, the simulated LHF was generally higher than buoy observations before 1200 UTC 15 September, especially at stations 2 and 5, but it then decreased quickly to the observed magnitude at station 4 after the appearance of evident SSC.

As two important factors determining heat flux, SST and surface winds are significantly affected by the surface status denoted by $ {C}_{\mathrm{d}} $ and $ {C}_{\mathrm{k}} $ (or $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $; Ming and Zhang, 2016). But the sensitivity of surface heat flux to these two parameters is different. Compared with the near-surface wind fields, the heat exchange was found to be more sensitive to $ {C}_{\mathrm{d}} $-induced SSC. With intense SSC, LHF was significantly inhibited around the core region of the typhoon, contributing to a weaker TC. However, when SST differed inconspicuously among the coupled simulations, the heat fluxes tended to be positively associated with the change of 10-m wind speeds. Therefore, different $ {C}_{\mathrm{d}} $ (or $ {C}_{\mathrm{h}} $/$ {C}_{\mathrm{d}} $) schemes won’t significantly affect the enthalpy fluxes through wind-induced processes, unless SSC is indistinctive. Wadler et al. (2021) and Zhang et al. (2015) also emphasized the importance of measuring atmospheric and oceanographic parameters in understanding TC intensity and structure changes. Moreover, the $ {C}_{\mathrm{d}} $ schemes examined here also affected the distributions of heat flux through continuous cooling of the deeper ocean, even long after the passage of TC. A larger $ {C}_{\mathrm{d}} $ promoted more momentum to be transported into the water columns, resulting in the continuous upwelling of the deeper colder water to the upper layer and the cooling of the OML (Fig. A3). This process was more evident to the right of the typhoon, consistent with the typically stronger wind field to the right of TCs. Therefore, a suitable $ {C}_{\mathrm{d}} $ scheme is essential for accurate numerical simulation of TCs.

The above results also indicate the complexity of choosing $ {C}_{{\rm{d}}} $ schemes in numerical simulations, especially for a moderate typhoon, as the threshold value of $ {C}_{\mathrm{d}} $ when becoming uniform or decreasing is uncertain. Powell et al. (2003) detected 33 m s⁻¹ as the threshold value. Vickery et al. (2009), however, suggested that the threshold value is up to 60 m s⁻¹. Since the intensity of Kalmaegi (2014) was weaker than 60 m s⁻¹, observations under stronger TC conditions were not obtained. Therefore, more numerical simulations and observations on the threshold values of $ {C}_{\mathrm{d}} $ should be conducted to draw a concrete conclusion on the speed of wind corresponding to $ {C}_{\mathrm{d}} $ leveling off or decreasing.

Another interesting phenomenon is the different intensification rates of the typhoon between the experiments with (the CSs cases) and without (the US cases) ocean feedbacks from 0000 UTC to 1800 UTC 15 September. Similar to the results found in earlier studies using an atmosphere-only model (Wang and Xu, 2010; Zeng et al., 2010), different $ {C}_{\mathrm{d}} $s in the US experiments rarely changed the intensification rates of the TC (dashed lines in Fig. 1b) but did lead to different maximum intensities of the mature TC. However, in the air–sea coupled model, the intensification rates of Kalmaegi (2014) were different using the three $ {C}_{\mathrm{d}} $ parameterization schemes when the wind speed became 35 m s⁻¹ or stronger. Both VMAX and the intensification rate in CS2 are undoubtedly the largest because of the decreased $ {C}_{\mathrm{d}} $ in high-wind conditions, while VMAX and the intensification rate in CS0 were the smallest due to the increasing sea surface roughness throughout the simulation. The results found here are different from earlier findings when using single atmospheric models (Wang and Xu, 2010; Zeng et al., 2010). This could be attributed to the enhanced enthalpy fluxes with complex air–sea conditions in the boundary layer. Suggested in previous studies (Zhang et al., 2017; Wadler et al., 2018; Zhang and Rogers, 2019), the surface heat flux-induced boundary layer recovery, which has been proven to regulate the thermodynamic structure of the boundary layer, plays an important role in TC development, including enhancing both the inflow and convergence and resisting the weak and dry downdraft near the core region, which in turn leads to stronger and more symmetric deep convection in TCs. Therefore, the role of SSC-impacted heat flux should be considered in future numerical simulations.

As is suggested in earlier studies, wave breaking and sea spray could significantly change the wave-current fields and enthalpy fluxes at high wind speeds, which are closely related to $ {C}_{\mathrm{d}} $ (Andreas and Emanuel, 2001; Moon et al., 2004; Bao et al., 2011). Therefore, a better understanding of the oceanic effects on $ {C}_{\mathrm{d}} $ during the early stage of a typhoon is needed to develop a suitable $ {C}_{\mathrm{d}} $ parameterization scheme in the coupled model. Besides, the discrepancies between in situ observations and model output also highlight the urgent need for advanced air–sea coupled models and more observation facilities in the future.

Acknowledgements. The authors are grateful to two anonymous reviewers for their constructive comments on this work and recommending several beneficial references. The authors thank Prof. Dake CHEN for providing the in situ mooring data, and Prof. Weimin ZHANG and Jianfang FEI for their help in improving the manuscript. This study has been supported by the National Natural Science Foundation of China under Grant Nos. 41775053, 41976003, and 42192552, and the National Key Research and Development Program of China under Grant Nos. 2019YFC1510001 and 2019YFC1510102. Additional support has been provided by the National Program on Global Change and Air-Sea Interaction (GASI-IPOVAI-04). We acknowledge model support from http://www2.mmm.ucar.edu/wrf/users/ and http://imedea.uib-csic.es/users/toni/sbpom/, provided by NCAR’s Antoni JORDI and Dong-Ping WANG, respectively. The FNL data for the WRF model was downloaded from https://rda.ucar.edu/datasets/ds083.2/. The ocean topographic data, ETOPO1, was acquired from the one-arc-minute global relief model of Earth’s surface (https://maps.ngdc.noaa.gov/viewers/wcs-client/). The ocean field for sbPOM was from https://www.HYCOM.org/data/glbu0pt08/expt-91pt1.

APPENDIX