Effects of Surface Flux Parameterization on the Numerically Simulated Intensity and Structure of Typhoon Morakot (2009)

Tropical cyclones (TCs) obtain heat and moisture from the ocean and transfer momentum back to the ocean at the air-sea interface. Thus, surface fluxes of sensible and latent heat play a very important role in the development and maintenance of TCs. (Malkus and Riehl, 1960) found that strong surface fluxes can increase the equivalent potential temperature in the TC eyewall, which is linked to the decrease in minimum sea level pressure at the TC center (i.e., increase in TC intensity). Early theoretical studies found that the TC intensity is sensitive to the selection of exchange coefficients for air-sea momentum and heat transfers (Ooyama, 1969; Rosenthal, 1971; Emanuel, 1986). (Rotunno and Emanuel, 1987) showed that the increase in the exchange coefficient for enthalpy transfer (C _K) and the decrease in the drag coefficient (C _D) caused the increase in TC intensity in idealized numerical models. Based on comparisons of model predictions with observations for a number of hurricanes, (Emanuel, 1995) concluded that the TC intensity is most sensitive to the ratio of C _K over C _D, and this ratio must lie in the range 0.75-1.5 to achieve consistent model simulations with observations. Full-physical numerical model simulations of TCs also suggest that the TC intensity is sensitive to C _D and C _K (e.g., Braun and Tao, 2000; Nolan et al., 2009).

Early observational studies (e.g., Smith, 1980; Large and Pond, 1981; Geernaert et al., 1988) investigated C _D in weak to moderate wind conditions with 10 m wind speed (U₁₀)<25 m s^-1 in the open ocean. They found that C _D increases almost linearly with U₁₀. The relationship between C _D and U₁₀ has the form 1000C _D=0.063U₁₀+0.61, where both C _D and U₁₀ have been adjusted to neutral stability (Smith, 1980). However, results from the Coupled Boundary Layer and Air-Sea Transfer (CBLAST) experiment indicate that C _D does not increase with U₁₀ with no limitation (Black et al., 2007; French et al., 2007). (Donelan et al., 2004) estimated C _D in a wave tank and found that C _D reaches a maximum for U₁₀ at ∼33 m s^-1 and then levels off at higher wind speeds. They attributed this leveling-off of C _D to the increasing sheltering of winds behind steeper waves at higher wind speeds. The above-mentioned studies support the result of (Powell et al., 2003), who estimated C _D using the profile method by fitting hundreds of GPS dropsonde wind profiles collected in hurricanes to a logarithmic relationship with height. (Ming et al., 2012) showed that C _D estimated based on observational data in typhoons also behaves similarly as in the (Powell et al., 2003) study.

Due to instrumentation limitations associated with sea salt and/or sea spray, direct measurements of surface enthalpy fluxes have been very difficult over the ocean. The only available in situ observations over the past 30 years (DeCosmo et al., 1996; Drennan et al., 2007) show that the exchange coefficient for latent heat transfer (C _E) does not depend on the wind speed. Similarly, the bulk sensible heat exchange coefficient (C _H), or Stanton number, has also been found to be independent of wind speed. Note that C _K is generally assumed to be equal to C _E and C _H, which is supported by the CBLAST result (Zhang et al., 2008). Wave tank observations also support the result from the CBLAST experiment, that C _K is nearly constant with U₁₀ up to ∼40 m s^-1 (Haus et al., 2010). (Bell et al., 2012) also found that values of C _K estimated based on budget analyses of dropsonde and radar data in the boundary layer of the eyewall region were close to those reported by (Zhang et al., 2008).

Most previous numerical studies on the impact of surface exchange coefficients on TC intensity and structure focused on Atlantic hurricanes or used idealized simulations (e.g., Montgomery et al., 2010; Bryan, 2012; Bao et al., 2012; Green and Zhang, 2013). How sensitive the simulated intensity and structure in West Pacific typhoons are to C _D and C _K is not well documented. The objective of this study is to investigate the effects of C _D and C _K on the simulated intensity and structure of Typhoon Morakot (2009). Emphasis is placed on simulations of structure, such as boundary layer heights, the warm core, and precipitation. Typhoon Morakot (2009) possessed unusual structures, including a large eye and strong outer circulation; thus, this storm was atypical compared to most hurricanes in the Atlantic basin. This motivates us to evaluate the choice of C _D and C _K in the numerical simulations to reproduce the unusual structures of Typhoon Morakot (2009).

Typhoon Morakot (2009) generated from a tropical depression over the Philippine Sea on 2 August 2009. It gradually intensified to a tropical storm and was assigned the name Morakot on 3 August. After that, it moved westward towards Taiwan. Morakot (2009) intensified to typhoon intensity at 1800 UTC 5 August, with the location of its center at (23.1^°N and 129.3^°E) and a maximum wind speed of 33.4 m s^-1. Morakot (2009) continued to move westward and made landfall in central Taiwan at 1200 UTC 7 August. One day later, Morakot (2009) turned to move northwestward and went back over water into the Taiwan Strait. Subsequently, Morakot (2009) weakened to a severe tropical storm and made its second landfall in Xiapu, Fujian Province, on the east coast of China at 0820 UTC 9 August (Fig. 1a). It gradually weakened as it moved northward. Roughly 24 hours later, it weakened to a depression during the landfall in Zhejiang Province, China. Typhoon Morakot (2009) brought record-breaking torrential rainfall over Taiwan (Chien and Kuo, 2011). This catastrophic storm sadly claimed more than 600 lives and caused more than 200 people to be classed as missing, creating a total cost of damage estimated at US $ 3.3 billion (Zhang et al., 2010). It also affected more than 11 million people throughout eastern China and damaged thousands of homes. As a result, Typhoon Morakot (2009) has received attention in several numerical studies.

Figure 1.  The (a) tracks and time series of simulated (b) maximum surface wind speed and (c) minimum central pressure (hPa) of Morakot (2009) from 0000 UTC 5 August to 0000 UTC 10 August 2009 from the three experiments and the corresponding best-track analysis by JMA.

(Schwartz et al., 2012) studied Typhoon Morakot (2009) using the Advanced Research Weather Research and Forecasting (WRF-ARW) modeling system with data assimilation of the microwave radiances using a cyclic, limited-area ensemble adjustment Kalman filter. They found that the track, intensity, and precipitation forecasts were improved after assimilating microwave radiances. (Xie and Zhang, 2012) performed ensemble simulations using WRF-ARW to investigate the dynamics and predictability of the record-breaking rainfall and flooding event in Taiwan induced by Typhoon Morakot (2009). They found that the large typhoon circulation and the southwesterly monsoon flow transported abundant moisture into southern Taiwan, which, along with the influence of the complex high terrain, produced the heavy rainfall.

(Wang et al., 2012) used the Cloud-Resolving Storm Simulator (CReSS) to study the dynamics related to the motion of Typhoon Morakot (2009) over Taiwan. Their simulations showed that the reduced moisture content induced the decrease in the rainfall and the increase in the storm translation speed. They also pointed out that the asymmetric precipitation in Typhoon Morakot (2009) played an important role in its very slow motion upon leaving Taiwan, the lengthening of the heavy-rainfall period, and the increase of the total rainfall amount. Furthermore, their results emphasized the potential contribution of asymmetric heating to the slowdown of the typhoon motion in the presence of complex terrain or in a monsoon environment. (Wang et al., 2013) demonstrated that CReSS successfully simulated and reproduced both the distribution and timing of the heavy rainfall in Taiwan with high accuracy when Typhoon Morakot (2009) passed by. The real-time forecast integrations of the CReSS model also showed high-quality quantitative precipitation forecasts.

(Hall et al., 2013) used the Advanced Regional Prediction System to simulate Typhoon Morakot (2009) when it made landfall over Taiwan. They investigated the mesoscale structure of Morakot (2009) and emphasized the role of deep convection on the rainfall simulation. They identified that relatively large-amplitude wave structures developed in the outer eyewall, known as vortex Rossby waves (VRWs). They found that the strong asymmetry of the convection was associated with wavenumber1 (WN1) VRWs, while the WN2 and WN3 VRWs were associated with the development of the deep convective band in Morakot's (2009) southwestern quadrant. (Huang et al., 2014) used the WRF model to explicitly simulate Typhoon Morakot (2009) and found that the simulated rain rate and precipitation efficiency (PE) over the Central Mountain Range (CMR) were highly correlated. They found also that the PE and the processes of vapor condensation and raindrop evaporation were strongly influenced by orographic lift. They also found that the increase in PE over the CMR compared with that over the ocean was due to an increase in the ice-phase deposition ratio when the liquid-phase condensation reduced as the air on the lee side subsided and moved downstream. They emphasized the effects of the terrain on the simulations of precipitation.

As mentioned above, the objective of this paper is to investigate the sensitivity of the simulated track, intensity and structure of Typhoon Morakot (2009) to the surface exchange coefficients. We used the ARW modeling system (version 3.2; Skamarock et al., 2008) to conduct numerical experiments. Three experiments were performed using the ARW model with two model domains. The model horizontal grid resolution was 4.5 km for the parent domain (D1) and 1.5 km for the nested domain (D2). The two domains covered 600× 400 and 421× 421 grid points, respectively (Fig. 1a). D2 was an automatic vortex-following moving nest and the center of the domain was always located at the center of the storm. In all the experiments, 35 vertical (σ) levels were used from the surface to the model top at 50 hPa.

D1 was initialized at 0000 UTC 4 August 2009 and was run for 6 days. D2 was started at 0000 UTC 5 August after a one-day spin-up. The outer domain was run in a four-dimensional data assimilation (FDDA) mode to provide the best possible large-scale conditions for the inner domain. The inner nested domain was run without FDDA. The Japan Meteorological Agency (JMA) 6-hourly gridded regional analyses at 20× 20 km horizontal resolution with 20 pressure levels were used as the initial and boundary conditions for the ARW model. The JMA analyses were produced using a multivariate three-dimensional optimum interpolation method to combine the first-guess fields from JMA's regional spectral model (RSM) with observations from a variety of platforms (JMA, 2007; Hosomi, 2005). Note that the JMA analyses did not capture the real structure and intensity of Morakot (2009). Compared to the best track, the error of the minimum sea-level pressure was 13.7 hPa, and that of the maximum surface wind speed was 10 m s^-1. The track of the storm in the JMA analyses also has a northwestward bias. Following (Ming et al., 2009), we appended a good vortex with the appropriate intensity (close to the best track) into the initial conditions. We ran D01 only from 0000 UTC to 1600 UTC 4 August, to make the storm intensity similar to the best track. The three-dimensional vortex was extracted and inserted back into the initial condition using the position of the best track, which corrected the bias of the storm center at the initial time. Under the new initial condition, the error of the minimum sea level pressure was 0.6 hPa and that of the maximum surface wind speed was 5.6 m s^-1.

The physics schemes used in the numerical experiments included the Purdue Lin microphysics scheme (Lin et al., 1983; Chen and Sun, 2002), the Yonsei University boundary layer scheme (Noh et al., 2003; Hong et al., 2006), the Rapid Radiative Transfer Model longwave radiation scheme (Mlawer et al., 1997), and the Dudhia shortwave radiation scheme (Dudhia, 1989). Cumulus parameterization schemes were not used for both domains.

In this study, three experiments (see Table 1) were designed to examine the impact of C _D and C _K on the typhoon intensity and structure simulations. As mentioned earlier, surface fluxes were parameterized through C _D and C _K. In the surface-layer schemes of WRF-ARW, C _D and C _K were parameterized using surface momentum (z₀) and thermal (z_{0 q}) roughness lengths. The drag coefficient in the neutral condition was defined as:

\begin{equation} \label{eq1} C_{ D}=\left(\dfrac{k}{\ln\frac{10}{z_0}}\right)^2 , (1)\end{equation}

where k is the von Kármán constant. The enthalpy exchange coefficient was defined as:

\begin{equation} \label{eq2} C_{ K}=\left(\dfrac{k}{\ln\frac{10}{z_0}}\right)\left(\dfrac{k}{\ln\frac{10}{z_{ 0q}}}\right) . (2)\end{equation}

In the control experiment (referred to as CTL hereafter), z₀ depended on the wind speed using the (Charnock, 1955) relationship, in the form of:

\begin{equation} \label{eq3} z_0=0.0185(u_\ast^2/g)+1.59\times 10^{-5} ,(3) \end{equation}

where u_* is the friction velocity, g is the gravitational acceleration. This formula causes z₀ to increase continuously with the wind speed, with no limit. For z_{0 q} the Carlson-Boland scheme (Carlson and Boland, 1978) was used:

\begin{equation} \label{eq4} \ln\dfrac{10}{z_{ 0q}}=\ln\left(\dfrac{10u_\ast k}{x_{ ka}}+\dfrac{10}{z_0}\right) , (4)\end{equation} where x _ka is a constant that equals 2.4× 10^-5. Since z_0q varies more slowly with the friction velocity than z₀, C _K increases more slowly with the wind speed than C _D (Dudhia et al., 2008). In the second experiment (referred to as TC1 hereafter), an alternate C _D formulation based on the high-wind laboratory experiment of (Donelan et al., 2004) was adopted. The roughness length was defined as: \begin{equation} \label{eq5} z_0=10\exp\left(-\dfrac{10}{u_\ast^{-1/3}}\right) , (5)\end{equation}

where 1.27× 10^-7≤ z₀ ≤ 2.85× 10^-3 (units: m). Values of C _D based on Eqs. (2) and (6) are smaller than those from the Charnok relationship. C _D increases almost linearly with the 10 m wind speed, up to a maximum of 0.0024 at 35 m s-1 (Davis et al., 2008). In addition, the thermal roughness length was calculated using the ramped formula (Dudhia et al., 2008) as follows:

\begin{equation} \ln\dfrac{10}{z_{0{ q}}}=\left\{ \begin{array}{l@{\quad}l} \ln(10^5) & u_{\ast}\le 1\\[2mm] \ln(10/[10^{-4}+10^{-3}(u_\ast-1)^2] & u_{\ast}\ge 1 \end{array} \right. . (6)\end{equation}

This setup of z₀ and z_{0 q} causes C _K to increase almost linearly with wind speed, and C _K becomes larger than C _D at wind speeds of >50 m s^-1. Note that the ratio C _K/C _D exceeds 1 as the wind speed becomes larger than 50 m s^-1, although the condition for the maximum surface wind speed of >50 m s^-1 was never met in our simulations (Fig. 1b). In the third experiment (referred to as TC2 hereafter), the drag formulation was still based on the result of (Donelan et al., 2004), but the thermal roughness length was adopted based on the Garratt formula (Garratt, 1992, p. 102) as follows:

\begin{eqnarray} \label{eq6} \ln\dfrac{10}{z_{0{ q}}}&=&\ln\dfrac{10}{z_0}+\ln\dfrac{z_0}{z_{0{ q}}} (7),\\ \label{eq7} \ln\dfrac{z_0}{z_{0{ q}}}&=&2.28R_{e^\ast}^{\frac{1}{4}}-2 ,(8)\\ \label{eq8} R_{e^\ast}&=&\dfrac{u_\ast z_0}{\gamma}(9) ,\\ \label{eq9} \gamma&=&[1.32+0.009(T-273.15)]10^{-5} , (10)\end{eqnarray}

where T is air temperature at the surface, γ is the kinematic viscosity of air (units: m² s^-1). In this setup, C _K is related to u_* and T. C _K decreases with increasing wind speed in smooth flow, while it is a constant in rough flow or high-wind conditions (Garratt, 1992, p. 102). Note that C _K using the Garratt formula is closer to recent observations from CBLAST and the wave tank experiment mentioned in the introduction. All simulations in the three experiments were run from 0000 UTC 5 August to 0000 UTC 10 August 2009, which covered almost the entire lifecycle of Typhoon Morakot (2009), from the intensifying stage to the landfalling stage. All experiments were conducted with the same initial vortex and boundary conditions, so the differences in the experiments were solely related to the different formulas of the exchange coefficients used in the WRF-ARW model.

The exchange coefficients and ratio of C _K/C _D as a function of the surface wind speed are shown in Fig. 2. In CTL, the Charnok relationship was used such that C _D increased almost linearly with the 10 m wind speed, because of the increasing surface roughness. C _K also increased slowly with the wind speed. In TC1 and TC2, when the Charnok relationship was changed to the new formulation of C _D based on (Donelan et al., 2004), C _D increased with the wind speed up to ∼33 m s^-1, then levelled off. In CTL, C _K was based on the Carlson-Boland formula and was larger than that based on the ramped formula used in TC1 for wind speeds of <30 m s^-1. On the other hand, the C _K quantities produced by these two formulas were close to one another for wind speeds of >30 m s^-1. The Garratt formula for the thermal roughness length was used in TC2, where the C _K was smaller than that in CTL and TC1 for wind speeds of >20 m s^-1. In TC2, the ratio of C _K/C _D was larger than that in CTL for wind speeds of <20 m s^-1, and marginally smaller for wind speeds of >20 m s^-1. Note that the ratio of C _K/C _D in TC1 was larger than that in the other two experiments for all wind speeds.

Figure 2.  The simulated (a) drag coefficient (C _D), (b) enthalpy exchange coefficient (C _K) and (c) exchange coefficient ratio (C _K/C _D) as a function of 10 m wind speed in the whole of domain 2 from the three experiments at 1800 UTC 6 August 2009.

The simulated track of Typhoon Morakot (2009) from the three experiments from 0000 UTC 5 August to 0000 UTC 10 August are compared with the JMA best track in Fig. 1a. As mentioned earlier, Typhoon Morakot (2009) moved westward in the first two days, and turned northwestward after making landfall in Taiwan. It then continued to move northward, and made a second landfall in Fujian Province. The simulated storm tracks in the three experiments were quite similar, although they deviated from the best track. All three experiments reproduced the observed west-northwestward storm motion. Overall, the simulated tracks are not sensitive to the different formulas of exchange coefficients.

The simulated maximum surface wind and minimum sea level pressure from the three experiments are compared with the best track from JMA in Figs. 1b and 1c. Although all three experiments were initialized with the same initial condition, the intensity forecasts were different, especially after the maximum wind speed reached 30 m s^-1 at the forecast time of 0000 UTC 6 August. In the next two days, all three experiments simulated maximum surface wind speed that was larger than observed. The maximum surface wind speed simulated in CTL and TC2 were closer to the best track. The difference in the simulated minimum sea level pressure, on the other hand, was clearer than that in the maximum surface wind speed among the three experiments. The simulated minimum sea level pressure in TC2 was closest to the best track among the three experiments.

Figure 3 plots the relationship between the maximum wind speed and minimum sea level pressure from the three experiments and the best track of JMA. It appears that the pressure-wind relationship is sensitive to C _D and C _K. The result is consistent with a recent numerical study carried out by (Green and Zhang, 2013), who found that the pressure-wind relationship was sensitive to the selection of different surface layer schemes in WRF simulations of Hurricane Katrina (2005). Our result is also consistent with (Bao et al., 2012), who ran different surface-layer physics in idealized simulations of the Hurricane Weather and Research Forecasting (HWRF) model and confirmed that the pressure-wind relationship was sensitive to C _D and C _K.

Figure 3.  Time-radius Hovmöller plots of azimuthally averaged tangential wind (units: m s^-1) at 2 km altitude from (a) CTL, (b) TC1 and (c) TC2, and radial wind (units: m s^-1) at 250 m altitude from (d) CTL, (e) TC1 and (f) TC2. The thick lines depict the RMW at 2 km altitude and the maximum inflow at 0.25 km altitude.

The simulated storm intensity and pressure-wind relationship in TC2 were closest to the best track among the three experiments, especially for wind speeds of >30 m s^-1. This result is encouraging because the C _D and C _K used in TC2 were closer to recent observations than those in the other two experiments. Our result generally indicates that when C _K is larger, the storm obtains more sensible and latent fluxes from the underlying ocean, resulting in larger surface wind speed and lower central pressure. On the other hand, when C _D is larger, surface friction is larger, which tends to reduce the surface wind speed. It was pointed out by (Montgomery et al., 2010) that larger surface friction can also lead to a larger gradient wind imbalance in the boundary layer, which could lower the central pressure. This explains why the simulated minimum sea level pressure in CTL was lower than that in TC2 but the maximum wind speed in CTL was smaller than that in TC2.

Figure 4.  Scatter plot of minimum sea level pressure vs maximum surface wind speed. The blue diamonds represent CTL, the red circles represent TC1, and the green crosses represent TC2. The open squares represent the pressure-wind relationship from JMA.

Time-radius Hovmöller diagrams of azimuthally averaged tangential wind velocity (V _t) at the altitude of 2 km and radial wind velocity (V _r) at the altitude of 250 m, from 0000 UTC 6 August to 1200 UTC 7 August, are displayed in Figs. 4a-c and 4d-f, respectively. As can be seen, the evolution of the axisymmetric V _t and V _r in the three experiments was similar. In the first 6 h of the simulations, Morakot (2009) was a weak storm that had a maximum V _t of 25 m s^-1 and minimum V _r of -10 m s^-1, located at the radius of ∼170 km. Later, V _t increased gradually with time and the radius of maximum tangential wind speed (RMW) became smaller. The magnitude of V _r doubled from -10 to -20 m s^-1 in the next 24 h and the radius of the peak inflow also contracted (Figs. 4d-f) in response to the intensification of the storm. Inward from the RMW, the radial velocity decelerated rapidly at a larger rate than outside the RMW. After another 16 h, the storm center was close to Taiwan such that the magnitudes of both V _t and V _r decreased. Among the three experiments, the peak V _t and V _r in TC1 were the largest, mainly because the ratio of C _K/C _D used in TC1 was larger than that in CTL and TC2.

The time-height evolution of the mean temperature anomalies (referred to as warm-core anomalies hereafter) from 0000 UTC 6 August to 1200 UTC 7 August is plotted in Fig. 5 for CTL, TC1 and TC2. Following (Liu et al., 1999) and (Li et al., 2013), we first calculated the mean temperature within the region of 630× 630 km from the typhoon's minimum surface pressure center at each vertical level. Then, the temperature anomaly was obtained by subtracting the mean temperature from the temperature at each grid point and each level. The warm-core anomaly was then defined as the average value of the temperature anomalies within the region of 300× 300 km from the storm center at each level. The height of the peak mean temperature anomaly represents the warm core height. It is evident from Fig. 5 that the peak warm core anomaly in CTL and TC2 was smaller than that in TC1, indicating that the warm-core anomaly was correlated with the storm intensity. According to the hydrostatic balance, the lower the minimum sea level pressure the larger the warm-core anomaly (Zhang and Chen, 2012). Thus, the warm-core anomaly in TC2 was the smallest among the three experiments as the storm intensity in TC2 was the lowest. On the other hand, the warm-core height was not correlated with the storm intensity, because the warm-core height was located at 8-12 km for the three experiments.

The horizontal distributions of 10 m wind speed valid at 1800 UTC 6 August 2009 are shown in Figs. 6a-c for CTL, TC1 and TC2, respectively. It is evident that the maximum surface wind speed in TC1 was larger than that in the other two experiments. The eyewall region, which covered the maximum wind speed, was broader in TC1 than in CTL and TC2. Comparing TC1 and TC2, the result suggests that increasing C _K alone would increase the storm intensity in terms of the maximum wind speed in the eyewall region. The surface wind distribution in TC2 was similar to that in CTL, but the maximum wind speed in the right-rear quadrant was slightly larger in TC2 than in CTL. Although C _K in TC1 was close to that in CTL, the maximum surface wind speed in TC1 was larger than that in CTL, especially for wind speeds of >25 m s^-1. This difference was mainly due to the fact that the C _D used in TC1 was smaller than that in CTL.

Large values of latent fluxes were found in the eyewall and primary rain band regions where surface wind speeds were also large (Figs. 6d-f). The maximum latent heat flux in TC2 was much smaller than that in the other two experiments, mainly because the C _K used in TC2 was smaller than that in CTL and TC1 (Fig. 2). The difference in the sensible heat fluxes among the three experiments (Figs. 6g-i) was much smaller than the difference in the latent heat fluxes, although the same formula was used for calculating both the latent and sensible heat fluxes. It appears that the simulated maximum sensible heat flux in TC1 was the largest among the three experiments, while that in TC2 was the smallest. The asymmetric distribution of the sensible heat flux in TC2 was closer to that in CTL than in TC1. Overall, the result (Fig. 6) implies that C _K and C _D alone have opposite effects on surface enthalpy flux, and C _K influences the enthalpy flux more than C _D. A larger C _K induces more sensible and latent fluxes, which support more energy for a storm to intensify. On the other hand, a larger C _D induces larger surface friction, which reduces the surface wind speed and in turn reduces the sensible and latent fluxes because these fluxes are also a function of the wind speed.

Figure 5.  Time-height diagrams of temperature deviation (units: K) from the three experiments: (a) CTL; (b) TC1; (c) TC2. The average was computed within the area of 300× 300 km from the surface minimum pressure center for simulations. The anomalies were obtained by subtracting the averaged temperature within the region at every height level.

Figure 6.  The model-simulated 10 m wind speed (color scale; units: m s^-1) and wind vectors (arrows; units: m s^-1) from (a) CTL, (b) TC1 and (c) TC2; latent heat flux (color scale; units: W m^-2) and wind vectors (arrows; units: m s^-1) from (d) CTL, (e) TC1 and (f) TC2; and sensible heat flux (color scale; units: W m^-2) and wind vectors (arrows; units: m s^-1) from (g) CTL, (h) TC1 and (i) TC2 at 1800 UTC 6 August 2009. The large vector indicates the motion of the storm, and the thin crosses divide the storm into four quadrants.

Figure 7.  Azimuthally averaged radius-height cross sections of tangential wind (units: m s^-1) from (a) CTL, (b) TC1 and (c) TC2; radial wind (units: m s^-1) from (d) CTL, (e) TC1 and (f) TC2; and the gradient force imbalance (F _PEC, units: m s^-1 h^-1) from (g) CTL, (h) TC1 and (i) TC2 at 1800 UTC 6 August 2009. The thick lines in (a-c) depict the height of the maximum wind speed varying with radius, and the thick dashed lines in (d-f) depict the inflow layer height, defined as the height where the radial wind speed is 10% of the peak inflow.

Figure 7 shows the radius-height plots of the azimuthally averaged tangential and radial velocities for all three experiments at 1800 UTC 6 August 2009. The magnitudes of the tangential and radial velocities in TC1 were generally larger than those in the other two experiments, consistent with the simulated storm intensities. However, the peak tangential wind speed in TC2 was smaller than that in CTL, which was not consistent with the intensity difference between these two experiments in terms of the maximum surface wind, as TC2 had a larger maximum surface wind speed. This discrepancy can be explained by the role of C _D in regulating the boundary layer dynamics. (Montgomery et al., 2010) pointed out that an increase in C _D leads to an increase in storm intensity in terms of maximum tangential wind speed in the boundary layer, although it reduces the surface wind speed through surface friction. Our result is consistent with that of (Montgomery et al., 2010), indicating that surface flux parameterization affects the vertical structure of wind velocities above the surface layer.

In many PBL schemes used in full-physics numerical models, one of the crucial elements is the boundary layer height, because it is coupled with the energy transport from the surface layer to the boundary layer and above (e.g., Beljaars and Viterbo, 1998; Noh et al., 2003). The boundary layer height is also a key variable that regulates the vertical distribution of turbulent fluxes and helps determine where turbulent fluxes tend to become negligible (Stull, 1988). Following (Zhang et al., 2011), the kinematic boundary layer height is defined by the height of maximum tangential wind speed (h _vtm). Inflow layer depth (\(h_\inf\)), defined as the height where the inflow reduces to 10% of the peak value, also represents the kinematic boundary layer height. In all three experiments, h _vtm decreased with decreasing radius toward the storm center (Figs. 7a-c). This behavior is consistent with the result given by (Zhang et al., 2011), who composited hundreds of dropsondes data collected from 13 hurricanes to study the characteristics of hurricane boundary layer heights. Within the radius of 125 km, h _vtm in TC2 was smaller than that in the other two experiments, and was closer to observations (Zhang et al., 2011, Fig. 5a).

Figure 8.  (a) U.S. Air Force Defense Meteorological Satellite Program (DMSP) polar-orbiting satellite F-17 microwave imagery of polarization corrected temperature with horizontal polarization with 91 GHz at 2112 UTC 6 August, and simulated composite radar reflectivity (units: dBZ) at 2100 UTC 6 August 2009 from (b) CTL, (c) TC1 and (d) TC2. (e) The observed composite reflectivity and simulated composite reflectivity from experiment (f) CTL, (g) TC1 and (h) TC2 at 1200 UTC 7 August 2009. (i) The observed composite reflectivity and simulated composite reflectivity from experiment (j) CTL, (k) TC1 and (l) TC2 at 0000 UTC 8 August 2009.

The difference in the inflow layer depth (\(h_\inf\)) among the three experiments was much larger than that in h _vtm (Figs. 7d-f). It appears that \(h_\inf\) in TC1 was the highest. All three experiments captured the decrease of \(h_\inf\) with deceasing radius, consistent with observations. Furthermore, h _vtm was smaller than \(h_\inf\) in all three experiments, which was also consistent with observations. Overall, the magnitude of \(h_\inf\) in TC2 was closest to observations (Zhang et al., 2011, Fig. 5b). The above results indicate that boundary layer heights are tied to the surface flux parameterization.

According to (Zhang et al., 2001), the momentum equation of the radial wind velocity in the cylindrical coordinates system can be written as

\begin{equation} \label{eq10} \dfrac{dV_r}{dt}=-\dfrac{1}{\rho}\dfrac{\partial p}{\partial r}+\dfrac{V_t^2}{r}+fV_t+2\Omega\cos\phi w\cos\lambda+U_{ d} , (11)\end{equation}

where w is vertical wind velocity; Ω is the angular velocity and φ is the latitude; r is the radius from the center; Λ is the azimuthal angle. Equation 11 states that the radial acceleration is determined by the radial pressure gradient force (F _P; first term on the right-hand side of the equation), the centrifugal force (F _E; second term), the Coriolis force (F _C; third and fourth terms), and diffusion (U _d; last term). The degree of gradient force imbalance or net agradient force (F _PEC) is evaluated by adding F _P, F _E, and F _C together (Figs. 7g-i). Firstly, F _PEC was larger in TC1 than in the other two experiments, supporting the fact that the simulated storm in TC1 was stronger, because the storm tended to spin-up faster when F _PEC was larger (Smith et al., 2009). Although the C _K/C _D was alike in CTL and TC2, F _PEC was larger in CTL than that in TC2. This was mainly due to the different C _D used in those two experiments. Following the dynamical explanation of (Montgomery et al., 2010), the agradient tendencies near the surface caused the inflowing rings of boundary-layer air to converge farther inwards in the storm center before rising out of the boundary layer and ascending into the eyewall updraught, resulting in enhanced maximum tangential wind speed. Our result is consistent with this argument (Fig. 7). It also suggests that intense positive supergradient acceleration occurs in the vicinity of the maximum tangential wind speed and is associated with the outflow jet above the boundary layer.

Next, we investigate the simulated radar reflectivity and precipitation in the three experiments. Figure 8a shows the observed U.S. Air Force Defense Meteorological Satellite Program polar-orbiting satellite F-17 microwave imagery of polarization corrected temperature with a horizontal polarization at 91 GHz valid at 2112 UTC 6 August. The simulated composite radar reflectivity valid at 2100 UTC 6 August from CTL, TC1 and TC2 are shown in Figs. 8b-d, respectively. Due to the interaction of the typhoon circulation with the monsoon flow and vertical wind shear (Wang et al., 2012), the storm became asymmetric. All three experiments captured the asymmetric distribution of radar reflectivity and reproduced the unclosed eyewall in the southern part of the storm. The simulated reflectivity in TC2 was only slightly closer to observations than in the other two experiments, as it captured broader high reflectivity area in the southern eyewall. Otherwise, the overall rainfall structure was similar in all experiments.

The observed and simulated radar reflectivity composites at 1200 UTC 7 August and 0000 UTC 8 August are shown in Figs. 8e-h and Figs. 8i-l, respectively. At 1200 UTC 7 August, Typhoon Morakot (2009) was located on the east side of Taiwan before landfall. It is evident from the observation (Fig. 8e) that the strongest reflectivity was located on the southwest side of the storm. This asymmetric rainfall pattern was captured by all the experiments. The simulated reflectivity in the eyewall region was stronger in CTL and TC1 than in TC2 (Figs. 8f-h). The simulated reflectivity in TC2 was slightly closer to observations than in the other two experiments because the unclosed eyewall, with the principal rainband located in the southern part of the storm and the strong echo over Taiwan, were captured in TC2. At 0000 UTC 8 August, the observed precipitation pattern became more asymmetric than in earlier periods. The strong echoes were observed on the south side of Taiwan (Fig. 8i). In CTL and TC1, the simulated reflectivity was still stronger than that in TC2 (Figs. 8j-l). Again, TC2 performed slightly better than in the other two experiments.

Overall, the total precipitation was strongest in TC1 and weakest in CTL, especially after 1200 UTC 6 August (Fig. 9a). Note that the accumulation period of precipitation is 1 h. As TC1 simulated the strongest storm while CTL simulated the weakest storm, this result suggests that the total precipitation is correlated with the storm intensity. The 10 m domain-averaged divergences of moisture flux (note that negative values represent convergence of the moisture flux) are shown in Fig. 9b. It is evident that the convergence of moisture flux was correlated with the total precipitation. This result is not surprising, as the low-level moisture was the main source of the rainfall. Nonetheless, the result suggests that surface flux parameterizations have a substantial impact on precipitation simulations.

Although the simulated rainfall over the ocean was strongly tied to surface flux parameterization, interestingly, we found that the rainfall over land (i.e., Taiwan) was much less sensitive to the surface flux parameterization. Figure 10 compares the 12 h accumulated precipitation from CTL, TC1 and TC2 together with the objective analyses of rainfall measured by automatic weather stations over Taiwan, valid at 0600 UTC and 1800 UTC 7 August. Prior to the landfall of Morakot (2009) in Taiwan (1800 UTC 6 August to 0600 UTC 7 August), the observation (Fig. 10a) shows two regions of strong precipitation across the island: one on the north side of Taiwan and the other on the south side over high mountains. All three experiments simulated these two regions of strong rainfall, although the simulated rainfall was much stronger than the observed valued. In particular, all three simulations over-predicted the precipitation from central to southern Taiwan. However, the difference in the precipitation distribution in the three experiments was very small over the whole island.

Figure 9.  Averaged (a) accumulated rainfall and (b) 10 m divergence of moisture flux from experiment CTL, TC1 and TC2. The average was computed within the area of 300× 300 km from the surface minimum pressure center.

Figure 10.  12 h accumulated precipitation (units: mm) valid at 0600 UTC 7 August 2009 from (a) automatic weather station hourly observations, (b) CTL, (c) TC1 and (d) TC2 over Taiwan; and 12 h accumulated precipitation valid at 1800 UTC 7 August 2009 from (e) automatic weather station hourly observations, (f) CTL, (g) TC1 and (h) TC2 over Taiwan.

The above result suggests that the precipitation was less sensitive to the exchange coefficients over land than over the ocean in Morakot (2009). Over the ocean, the role of the exchange coefficient in precipitation was an indirect one, through influencing the (horizontal) moisture flux convergence (and ultimately rainfall) and through affecting the intensity of the storm. For the precipitation over land, the moisture came from the ocean, even though the end precipitation fell over land. In Typhoon Morakot (2009), the effect of terrain (i.e., forced uplift) played a dominant role in the distribution of the precipitation over Taiwan (Hall et al., 2013; Wang et al., 2013), which is likely the main reason for the similar rainfall simulations among the three experiments.

In this study, three numerical experiments were performed with the WRF-ARW model to study the impact of surface flux parameterizations on the structure and intensity of Typhoon Morakot (2009). The initial conditions of the three experiments were all from the JMA RSM analysis field. The simulated track and intensity of Morakot (2009) were verified against the best track. Different formulas of momentum and heat roughness lengths were tested in sensitivity experiments that governed the behavior of the surface exchange coefficients for momentum and heat transfers. The results showed that the simulated track was not sensitive to the exchange coefficients, but the simulated intensity and structure were.

Our results indicate that the surface exchange coefficients are key factors for the simulation of surface wind speed and fluxes. The effect of C _K on the surface enthalpy flux is straightforward because of the linear relationship between these two parameters. On the other hand, the effect of C _D on the enthalpy flux takes place via the surface wind speed. When the C _D is small, the maximum surface wind speed tends to be larger due to reduced surface friction. In turn, the enthalpy flux becomes larger because of the larger wind speed. Overall, we found C _K had a larger impact on the enthalpy flux simulation than C _D.

Consistent with previous studies (e.g., Emanuel, 1995), the simulated storm intensity was found to be more sensitive to the ratio of C _K/C _D than to C _K or C _D alone. According to the idealized numerical simulation given by (Montgomery et al., 2010), C _K/C _D should have a critical value for the intensification of storm. If C _D is too large, the storm will not intensify. When C _K/C _D is larger, the simulated storm is stronger and vice versa. In the CTL experiment, the intensity simulation is comparable to that in TC2, because a similar C _K/C _D was used in these two experiments.

The pressure-wind relationship was also found to be sensitive to C _D and C _K, consistent with recent numerical studies of Atlantic hurricanes (Bao et al., 2012; Green and Zhang, 2013). Overall, the simulated intensity and pressure-wind relationship in TC2 was closest to the best track than those in CTL and TC1. This result is encouraging because the C _K and C _D used in TC2 were close to recent field and wave tank observations. This result is also consistent with that of (Zhang et al., 2012), who showed that observation-based surface layer and boundary layer physics led to improvements in the operational HWRF model and better intensity forecasts.

Our results also indicate that simulated structures, such as the surface wind distribution, boundary layer heights, warm-core anomaly and height, and precipitation are affected by C _D and C _K. Compared to the dropsonde observations from (Zhang et al., 2011), the simulated kinematic boundary layer heights in TC2 are closer to observations than the other two experiments. The warm-core anomaly is tied to the storm intensity but not the warm-core height, consistent with Stern and Nolan (2012). The difference in the rainfall over the ocean is consistent with the difference in storm intensity, which can be explained by the difference in the convergence of moisture flux in the boundary layer. Over land, the simulated rainfall is much less sensitive to C _D and C _K than over the ocean, which we attribute to the dominance of the terrain effect on the precipitation in Typhoon Morakot (2009), as pointed out by (Wang et al., 2013) and (Hall et al., 2013).

We also conducted dynamical analyses to investigate why C _D and C _K affect the vertical structure of wind velocities in the boundary layer. Consistent with (Montgomery et al., 2010), we found that a larger drag coefficient can lead to a larger gradient wind imbalance in the boundary layer (Fig. 7). As a result of the larger agradient forcing, the boundary-layer air converged farther inward near the storm center before rising out of the boundary layer and ascending into the eyewall updraft. The end result was enhanced maximum tangential wind speed, despite the loss of absolute angular momentum en route.

In this study, we focused on investigating the sensitivity of the simulated intensity and structure of Typhoon Morakot (2009) to the surface exchange coefficients only, while keeping the rest of the model physics the same. We note that other parts of the model physics (e.g., planetary boundary layer parameterization and radiation parameterization) may also be important for TC simulations. Future work will evaluate the impact of other aspects of model physics on numerical simulations of TC structure and intensity change.