In this study, the Weather Research and Forecasting Model (WRF)–Advanced Research WRF (ARW) modeling system (version 3.8.1; hereafter, WRF-ARW) was used (Skamarock et al., 2008). In all simulations, all physics options other than the surface flux options (described below) were identical, including the WRF single-moment three-class microphysics scheme (Hong et al., 2004), updated Kain–Fritsch scheme (Kain and Fritsch, 1990; Kain, 2004) with a moisture-advection-based trigger function (Ma and Tan, 2009), Rapid Radiative Transfer Model for Global Circulation Models (RRTMG) shortwave and longwave schemes (Iacono et al., 2008), Yonsei University boundary layer scheme (Hong et al., 2006 for YSU) used with the revised fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5) similarity surface layer scheme (Jiménez et al., 2012), and the unified Noah land surface model (Tewari et al., 2004) over land. In WRF, the bulk transfer coefficients as well as all surface fluxes (momentum, sensible heat, and latent heat) are parameterized in a surface layer scheme. A brief description of the surface flux parameterization in WRF is provided in Appendix A for reference.

2.1 Flux parameterization in WRF-ARW

The behaviors of bulk transfer coefficients are strongly dependent on the momentum and scalar (heat and moisture) roughness lengths (see Eqs. A7, A8, and A9). Over the ocean, roughness length depends on the varying surface sea state, which is considerably ascribed by wind waves (e.g., Geernaert et al., 1987; Smith, 1988; Fairall et al., 2003). Under high-wind conditions, wave breaking and wind tearing the breaking wave chest generate sea foam and sea spray, thereby leading to significant regime change in the surface sea state, which in turn affects the specification of the roughness length (e.g., Powell et al., 2003; Donelan et al., 2004; Zweers et al., 2010). For the surface layer physics over the ocean, WRF-ARW provides alternative formulations of aerodynamic roughness lengths of the surface momentum and scalar fields for tropical storm application, which can be set through the namelist variable isftcflx. There are three isftcflx options available in WRF-ARW. Here, we study the resolution dependence of all three. Each flux option is referred to by its corresponding number (0, 1, and 2) in the WRF namelist file.

Subsequently, we describe the three sets of formulas of roughness lengths for momentum, sensible heat, and latent heat that have been implemented in the revised MM5 scheme in WRF-ARW. For clarity, we calculated the bulk transfer coefficients and the momentum scaling parameter (friction velocity) with respect to the three sets of implemented roughness length formulas in the case of neutral stability and a reference height of 10 m, using a set of pseudo-wind speeds. Under neutral condition, the stability parameters $(z+z_{\mathrm{0}})/L$ and $(z+z_h)/L$ are zero (see Eqs. A7, A8, and A9); therefore, all terms that involving the stability function can be dropped. The results are depicted in Figs. 1 and 2.

Figure 1Plots of (a) roughness lengths and (b) bulk transfer coefficients versus wind speed at 10 m height for neutral conditions. Flux schemes F0, F1, and F2 are colored black and gray, red and pink, and green. The thick solid curves are z₀ and C_D, thin solid curves are z_h and C_H, and dashed curves are z_q and C_Q. It is noted that some curves are identical. For F1 and F2, z₀ (and thus C_D) is the same; for F1 and F2 the respective heat and moisture roughness and coefficients are the same: z_h=z_q and C_H=C_Q.

Download

Figure 2Plots of (a) momentum transfer coefficients, (b) friction velocity, (c) scalar transfer coefficients, and (d) transfer coefficient ratios versus wind speed at 10 m height from available observational estimates as found in the literature. The scalar transfer coefficients shown here can be C_H, C_Q, or C_K upon the availability in the literature that we referred to in this study. The purple dots represent the values from Bell et al. (2012), blue circles represent values from French et al. (2007), dark green circles represent values from Geernaert et al. (1987), green curve represents values from Jarosz et al. (2007), red curve represents values from Large and Pond (1981), red dots represent values from Powell et al. (2003), blue dots represent values from Richter et al. (2016), and the orange dots represent values from Vickers et al. (2013). The curves are given by fitted equations; the circles represent the available data points; the dots represent average values in a range of wind speed bins with the bars showing the corresponding data spread as denoted by either the standard deviations or the 95 % confidence intervals, where the wind speed bin and the representation of data spread vary in the literature. Table B1 provides brief information of these observational estimates. The respective bulk parameters derived from WRF, represented as shown in Fig. 1, are superimposed for a comparison.

Download

2.1.1 Momentum and moisture roughness length for flux option F0

For the default flux option (F0), the momentum roughness length is given as Charnock's (1955) expression plus a viscous term, following Smith (1988):

$$\begin{array}{}\text{(1)}& z_{\mathrm{0}}=\mathit{\alpha }{\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{g}}+{\displaystyle \frac{\mathrm{0.11}v}{u_{\ast }}},\end{array}$$

where α is the Charnock coefficient and v the kinematic viscosity of dry air, for which a constant value of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ is used. The Charnock expression relates the aerodynamic roughness length to friction velocity for describing the gross effect of wavy sea surface, which is generated and supported by wind stress; the viscous term describes roughness behavior under smooth flow conditions. The use of this roughness length formula yields a monotonic increase in the momentum transfer coefficient, as well as roughness length, with wind speed exceeding approximately 3.0 m s⁻¹ (Fig. 1; black thick solid curve). A constant value of α=0.0185 is used for F0.

The moisture roughness length is adapted from COARE 3.0 (Fairall et al., 2003) and expressed as a function of the roughness Reynolds number, with an additional lower limit of $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{9}}$ used in this option:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{z}_q=& \max \left[\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{9}},\right.\\ & \min \left.\left(\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}},\mathrm{5.5}\times {\mathrm{10}}^{-\mathrm{5}}{R}_{\ast }^{-\mathrm{0.6}}\right)\right],\end{array}\end{array}$$

where $R_{\ast }=z_{\mathrm{0}}{u}_{\ast }/v$ is the roughness Reynolds number and v the kinematic viscosity of dry air and a function of air temperature for the R_∗ calculation. The temperature roughness length z_h is set to equal to z_q. The upper limit of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ keeps the scalar roughness length invariant at extreme wind speeds (below approximately 4.0 m s⁻¹), whereas the lower limit of $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{9}}$ is implemented to prevent the model from blowing up. Below this high wind speed, the heat and moisture transfer coefficients increase monotonically with wind speed (Fig. 1; gray thin solid and black dashed curves). For this option, the formulas of roughness lengths are different from that used by Green and Zhang (2013), where the viscous term for z₀ was given as a constant value of $\mathrm{1.59}\times {\mathrm{10}}^{-\mathrm{5}}$, z_h is set equal to z₀, and z_q is set to a somewhat different functional form of z₀ and u_∗.

2.1.2 Momentum and moisture roughness length for flux option F1

For the flux option 1 (F1), the momentum roughness length is expressed as a blend of two roughness length formulas (Green and Zhang, 2013):

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}{z}_{\mathrm{0}}=& max\left[\mathrm{1.27}\times {\mathrm{10}}^{-\mathrm{7}},\right.\\ & min\left.\left(z_{\mathrm{w}}{z}_{\mathrm{2}}+\left(\mathrm{1}-z_{\mathrm{w}}\right)z_{\mathrm{1}},\mathrm{2.85}\times {\mathrm{10}}^{-\mathrm{3}}\right)\right],\end{array}\text{(3b)}& {\displaystyle }{z}_{\mathrm{w}}=min\left(\mathrm{1},{\left[{\displaystyle \frac{u_{\ast }}{\mathrm{1.06}}}\right]}^{\mathrm{0.3}}\right),\text{(3c)}& {\displaystyle }{z}_{\mathrm{1}}=\mathrm{0.011}{\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{g}}+\mathrm{1.59}\times {\mathrm{10}}^{-\mathrm{5}},\text{(3d)}& {\displaystyle }{z}_{\mathrm{2}}={\displaystyle \frac{\mathrm{10}}{\exp \left(\mathrm{9.5}{u}_{\ast }^{-\mathrm{1}/\mathrm{3}}\right)}}+{\displaystyle \frac{\mathrm{0.11}v}{max\left(u_{\ast },\mathrm{0.01}\right)}},\end{array}$$

where the first roughness length formula (z₁) is again the Charnock (1955) expression plus a constant value of the viscous term, with α=0.011, as suggested by Smith (1988). The second formula (z₂) is the exponential expression from Davis et al. (2008) plus a viscous term; here a constant kinematic viscosity $v=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ is used. The two roughness length formulas are combined using a weight function (z_w), with a lower and upper limit on z₀ of $\mathrm{1.27}\times {\mathrm{10}}^{-\mathrm{7}}$ m and $\mathrm{2.85}\times {\mathrm{10}}^{-\mathrm{3}}$ m, respectively. The lower and upper limits on z₀ are adapted from Davis et al. (2008), with a slightly different value of the lower limit. The upper limit used here acts to prevent a monotonic increase in the momentum roughness length as well as the momentum transfer coefficient, under high-wind conditions. The resulting curves present an increase with wind speed below 33.0 m s⁻¹ but exhibit a leveling-off behavior starting at the wind speed of 33.0 m s⁻¹ (Fig. 1; red thick solid curves). The wind speed of 33 m s⁻¹ is extremely close to the lower limit of maximum wind speed for a typhoon. The use of this z₀ formula results in apparently smaller momentum roughness length and transfer coefficient at all wind speeds than that in F0. This leveling-off in C_D with wind speed suggests a decline in the efficiency of the exchange of momentum across the air–sea interface (see Eq. A1). Considering that the drag effect acts as a momentum sink in the surface layer, it is expected that an atmospheric simulation using the z₀ formula of F1 shall retain a larger portion of the total amount of near-surface momentum at all wind speeds than that using the F0 formula, and such a difference in the transfer efficiency becomes increasingly significant with higher wind speeds because of the leveling-off behavior seen in F1.

F1 sets $z_h=z_q={\mathrm{10}}^{-\mathrm{4}}$ m for all wind speeds. Large and Pond (1982) suggested constant values for z_h and z_q under different conditions of stability, ranging from 10⁻⁹ m (stable) to 10⁻⁴ m (unstable). The constant value of 10⁻⁴ m used here corresponds roughly to the unstable condition in Large and Pond (1982). While using an invariant scalar roughness length, the heat and moisture transfer coefficients can still vary with wind speed (Fig. 1; pink thin solid and red dashed curves) because of the drag-dependent effect (see Eqs. A7, A8, and A9). Consequently, the heat and moisture transfer coefficients are also beginning to level off at the wind speed of 33.0 m s⁻¹.

2.1.3 Momentum and moisture roughness length for flux option F2

For the flux option 2 (F2), the momentum roughness length is the same as that for F1. The temperature and moisture roughness lengths are expressed based on the formula proposed by Brutsaert (1975a):

$$\begin{array}{}\text{(4)}& {\displaystyle }{z}_h=z_{\mathrm{0}}\exp \left[-\mathit{\kappa }\left(\mathrm{7.3}{\mathrm{R}}_{\ast }^{\mathrm{1}/\mathrm{4}}{P}_r^{\mathrm{1}/\mathrm{2}}-\mathrm{5}\right)\right],\text{(5)}& {\displaystyle }{z}_q=z_{\mathrm{0}}\exp \left[-\mathit{\kappa }\left(\mathrm{7.3}{\mathrm{R}}_{\ast }^{\mathrm{1}/\mathrm{4}}{S}_c^{\mathrm{1}/\mathrm{2}}-\mathrm{5}\right)\right],\end{array}$$

where R_∗ is the roughness Reynolds number, κ is the von Karman constant, P_r and S_c are the Prandtl and Schmidt numbers, and 7.3 and 5 are experimental constants (Brutsaert, 1975b). As in F0, a temperature-dependent v is also used for R_∗ calculation. The values of κ, P_r, and S_c are set to 0.40, 0.71, and 0.60, respectively, as suggested by Brutsaert (1979). This set of formulas relates the heat and moisture roughness length to a scaling of z₀ by exponential expression, where the exponent contains information about the near-surface flow regime as characterized by the roughness Reynolds number and Prandtl and Schmidt numbers. The resulting z_h and z_q are practically an exponential decay of z₀ with wind speed (Fig. 1; green thin solid and green dashed curves). The scalar roughness lengths in F0 and F2 both decrease with wind speed, and their decreasing rates as well as magnitudes are similar at higher wind speeds. However, their resulting scalar transfer coefficients are apparently different. Unlike the monotonic increasing behavior found in F0, for F2 the heat and moisture transfer coefficients show an upward trend with wind speeds below 33.0 m s⁻¹ and then start decreasing at this wind speed. Such a discrepancy between F2 and F0 results from the drag-dependent effect in the scalar transfer coefficients (see Eqs. A7, A8, and A9).

2.3 Experimental designs

Table 1 summarizes the experiments designed for this study. There were three resolution groups, with horizontal grid spacing of 1, 3, and 6 km. All three isftcflx surface flux options were tested for these three resolution groups. We considered 3kmF0 as the control experiment. All experiments used a single domain and 45 vertical levels. The model domain covered a spatial region bounded by 0–19.62^∘ N latitudes and 117.88–147.12^∘ E longitudes (Fig. 3; red box). The experiment procedure included a nudging and a sensitivity stage. For the nudging stage, the grid nudging (Staufer and Seaman, 1994) implemented in WRF was applied to each resolution group using the default flux option (F0). The numerical solutions were nudging toward the gridded reanalysis for 24 h starting at 00:00 UTC on 4 November 2013 (Fig. 3; tan dots). In this study, the analysis nudging was applied to the horizontal wind components, potential temperature, and water vapor mixing ratio. The nudging coefficients for all variables were set at 0.0003 s⁻¹. Nudging was only applied at all levels above the PBL to allow for the development of mesoscale processes near the surface. The National Centers for Environmental Prediction Global Forecast System (GFS; Environmental Modeling Center 2003, 2016) operational global analysis dataset was used to provide the initial and boundary conditions. The GFS global analysis has a horizontal resolution of 0.5^∘, and it is also used as data for the grid nudging. The model sea surface temperatures (SSTs) were taken from the RTG-SST dataset (Gemmill et al., 2007), with daily time resolution; there was no ocean feedback in this study. The SSTs varied during the nudging stage up until 00:00 UTC on 5 November 2013 and thereafter remained fixed throughout the sensitivity simulation such that the differences in underlying SSTs between the flux options became negligible. Consequently, sensitivity of typhoon intensity and structure to the flux options (thus the surface flux formulas) can be determined to a large extent. All the sensitivity experiments started at 00:00 UTC on 5 November 2013 and ended at 06:00 UTC on 8 November 2013 (Fig. 3; red dots). All results are from within the 78 h that the experiments ran for.

Figure 3The model domain (red box) and calculation domains for the contoured frequency by altitude diagrams (blue box) and for energy spectra (cyan box). The red, tan, and green and dots represent the locations of Haiyan obtained from the best-track dataset. The red dots indicate the period of sensitivity experiments, the tan dots indicate the nudging period, and the green dots are the locations outside our simulation period.

Download

Table 1Summary of the resolutions and surface flux options tested. The nudging period for all experiments is between 00:00 UTC on 4 November 2013 and 00:00 UTC on 5 November 2013, with an identical nudging setting for each resolution group using the default flux option (F0). Thereby the corresponding sensitivity experiments start with an identical initial TC (tropical cyclone) structure and location. The MCP and MSW indicate, respectively, the minimum central pressure at the TC center and maximum surface wind speed near the TC center.

Download Print Version | Download XLSX

Typhoon Haiyan formed in the western North Pacific on 3 November 2013 and underwent a rapid intensification a few days later. Haiyan reached peak intensity of 895 hPa before it made a cataclysmic landfall near Guiuan, Eastern Samar, Philippines, on 8 November 2013 at about 04:40 PHT (UTC + 8). The storm then crossed the Philippines, migrated west-northwestward, and made its final landfall in Vietnam. Our simulation period for the sensitivity experiments covers the stages of Haiyan's peak intensity and landfall. We focus on the mature stage of Haiyan before its landfall near Guiuan, Eastern Samar.

All simulations show a similar west-northwestward movement toward land and appropriately follow the observed best track, but with slower translation speeds (Fig. 4a). The typhoon landfall at Samar island is delayed between 3 and 5 h in these sensitivity simulations. Because of the delayed landfall, we compare the evolution of the typhoon related to its longitudinal position (Fig. 4b). In general, the simulations continue to intensify and reach their peak right after 18:00 UTC on 7 November 2013 and before landfall. We refer to the time span from 18:00 UTC on 7 November 2013 to 00:00 UTC on 8 November 2013 as the mature stage of the simulated typhoon. In this study, the cumulus parameterization was used for all resolutions. This is because we intend to keep consistency among all cases. Results from an additional simulation of the 1kmF0 with the cumulus parameterization turned off revealed that the simulated Haiyan intensity and structure are overall similar. The minimum central pressure and storm track of the 1kmF0 without cumulus parameterization follow almost the same evolution as that of 1kmF0 (Fig. S3). In our case, the use of a cumulus parameterization in high-resolution (1 km) simulation does not exert considerable influence on the typhoon intensity, track, and structure. Some studies have also revealed that the activation of cumulus parameterization for simulation with grid spacings of 2–3 km produced an overall similar simulated storm to the simulation with explicit convection (e.g., Yu et al. 2011; Li et al. 2018; On et al. 2018).

Figure 4(a) Haiyan's best track (black dots) and simulated tracks of 1 km (blue lines), 3 km (red lines), and 6 km (green lines) from 00:00 UTC on 5 November 2013 through 06:00 UTC on 8 November 2013. Thick and thin curves represent a span of 24 h alternately. (b) The evolutions of the best track and simulated minimum central pressures (hPa) shown in relation to longitude. The longitudinal position of the first landfall of Haiyan is around 125^∘ E. The gray circle on each curve denotes the corresponding data point at 18:00 UTC on 7 November 2013. (c) Scatterplot of maximum wind speed versus minimum central pressure. Reference curves Koba and AH_KZ are the operational WPRs for the JMA and JTWC, respectively. Koba_1min is the modified Koba curve with 10 min sustained winds converted to 1 min averaged values. The detailed tracks and central positions of the simulated storms during their mature stage are given in Figs. S1 and S2 in the Supplement.

Download

For each resolution group, F1 is the most intense in terms of minimum central pressure, whereas the intensity of F2 is somewhat between that of F0 and of F1. The intensity is sensitive to model resolution. Overall, the intensity increases as the resolution is changed from 6 to 3 km, but it significantly increases as grid spacing is reduced to 1 km. The experiment 3kmF1 yields an intensity level comparable with that of 1kmF0, demonstrating the benefit of using a more reasonable flux option in improving intensity simulation. Accordingly, the combination of the decrease in the momentum transfer coefficient (generally tends to reduce the energy loss) and increase in enthalpy transfer coefficients (thus more energy gain) has greater potential to yield a stronger storm. Finally, a higher resolution is more conducive to the intensification due to the change in flux formulas.

The comparison between the simulated maximum winds near the typhoon center and the best-track wind information is somewhat complicated owing to the wind speed averaging period (Kueh, 2012; Knapp et al., 2013). We use the International Best Track Archive for Climate Stewardship (IBTrACS) dataset for comparison, which compiles best-track information from various agencies worldwide (Knapp et al., 2010). The best-track information in Figs. 3 and 4a is taken from the WMO subset of the IBTrACS (IBTrACS-WMO, v03r09). The WMO best track for Haiyan is taken from the best-track data provided by the Japan Meteorological Agency (JMA), which relies on a wind pressure relationship (WPR) where the maximum wind is given in terms of 10 min sustained winds (Koba et al., 1991). However, simulated maximum winds near the typhoon center were derived from WRF hourly output instantaneous maximum 10 m wind speeds at the beginning of each hour. Figure 4c presents the scatterplot of maximum wind speed versus minimum central pressure. The discrepancy found between the simulations and the best-track data may be partially due to different wind speed averaging periods. The WMO best-track wind pressure distribution appropriately follows the operational WPR, namely the reference curve Koba. We use a factor of 0.88, which is used operationally at the Joint Typhoon Warning Center (JTWC), to convert the 10 min sustained wind in Koba's WPR to an equivalent 1 min sustained wind. The modified Koba curve Koba_1min now agrees with the operational WPR used at the JTWC. The JTWC's WPR uses 1 min sustained winds (Atkinson and Holliday, 1977; Knaff and Zehr, 2007). In practice, operational WPRs were derived statistically to best fit a set of observed parameters (pressures and winds) by using the gradient wind approximation as a fundamental guide. In the context of statistical regression and use of observed parameters, other information, such as the vortex size, ambient conditions, and geographical location, would be incorporated into these reference curves (Knaff and Zehr, 2007). The operational WPR may somewhat be considered the best fit to the actual mass–wind relation of TCs. The winds and pressures of Haiyan from the JTWC best-track data have been plotted in Fig. 4c, and the distribution agrees well with the JTWC operational WPR; we thereby omitted them from the plot. All our simulated wind pressure distributions present similar slopes to the two reference curves of 1 min averaged winds, indicating that these simulations are valid for further investigation of the typhoon structures.

Bao et al. (2012) reported that the simulated WPR is controlled by the formulation of the drag coefficient. From a set of simulations for North Atlantic TCs between 2008 and 2011, Green and Zhang (2013) found that different C_D formulas affect wind pressure distribution (see Fig. 9 in Green and Zhang, 2013). We found insufficient evidence supporting any apparent slope change in wind pressure distributions of different flux options. This discrepancy may be partially due to the use of the square of the maximum 10 m wind speed in their wind pressure distribution plot. By comparing the respective linear regression curves between all experiments, we found that the differences in flux options and resolution result in a slight change in the slope (data not shown). However, no distinct dependence of the slopes on the choice of flux option and resolution could be identified. In our study, the differences in the wind pressure distribution between these sensitivity simulations are demonstrated by their magnitudes. For each resolution group, F1 is the most intense in terms of both minimum central pressures and maximum winds, whereas the intensity of F2 is somewhat between that of F0 and F1. A comparison of resolution groups used in F0 reveals a significant increase in intensity at higher resolutions during the mature stage (where the storm has a higher strength), with 1kmF0 receiving larger increments. The set of experiments using F1 presents a similar pattern, with the larger increments in 1kmF1 being much more marked.

In brief, increased resolution and more reasonable flux parameterization (e.g., F1 in this case) can both improve typhoon intensity simulation to a certain extent. However, a sufficiently high resolution is more likely to the benefit from improved flux formulas. The implication is 2-fold. First, the use of reasonable flux formulas to improve intensity simulation is perhaps more efficient than using an extremely high resolution when computational resources are limited. Second, higher model resolution is conducive to improving surface flux representation in strong typhoon intensity simulations. The mechanisms underlying these two issues are notable. We address the issue of higher model resolution and its effect on surface flux representation in the following sections.

The sensitivity of simulated typhoon intensity to various resolutions and flux options is also associated with changes in the simulated typhoon structure. To isolate differences between simulations of different grid spacing, the hourly outputs on the WRF native grid were transformed to the cylindrical (polar-height) coordinate system centered at the simulated typhoon at each time point. To avoid spatial sampling issues in the comparison of mass and wind fields, all model outputs were transformed to an identical resolution setting. We used a 3 km horizontal resolution and an uneven vertical resolution with smaller spacing in the boundary layer for the cylindrical transformation. Specifically, the cylindrical coordinate used in this study has an azimuthal resolution of 2^∘, a radial resolution of 3 km, and an uneven vertical resolution within the radius and height of 900 and 18 km, respectively. All the data to be analyzed are interpolated from the model native coordinate to the cylindrical coordinate using bilinear interpolation. The zonal and meridional components of horizontal wind fields are transformed to tangential and radial wind components using their standard formulas. The (height-invariant) typhoon center is defined as minimum central pressure using the surface pressure field. We also examined the geometric center for the local minimum pressure and found that the distances between the center at grid point and the geometric center are less that the half of the grid spacing, for the simulation period after 24 h (valid at 00:00 UTC on 6 November 2013). The magnitudes of the vertical tilt of the center (in terms of the WRF output perturbation pressure field) are also small: about the size of one grid box. In this study, radial variations of model fields were all derived from these transformed subsets.

Figure 5 shows the horizontal distributions of wind speed at 10 m height, total rainfall, and surface enthalpy fluxes for F0 with 3 km resolution after the cylindrical transformation. The wind speed distribution shows a ring of higher wind speed (exceeding approximately 42 m s⁻¹), with wider regions of high wind speed over the right side (northern semicircle) of the track, with a partial ring of the maximum wind speed (∼50 m s⁻¹), open to the west, at the range of radii 40–50 km. A region of intense rainfall is located on the left-rear side of the track. Regions of a larger amount of surface fluxes are located near the local maximum wind speeds; however, they are not necessarily co-located. This is because the near-surface enthalpy vertical gradient, other than wind speed, also influences flux calculation. The discrepancies in the horizontal distributions of surface fields between different resolution groups are more significant than those between flux options (Fig. 6). The 1 km resolution group yields a significant eyewall (in terms of these surface fields) contraction compared with the other two groups. The 3 km resolution group exhibits a relatively larger eye region, less contracted radius of the maximum wind speed, and total rainfall and surface fluxes. The 6 km resolution group does not present a well-defined eyewall structure in terms of these surface fields. Therefore, the increase in resolution from 6 to 3 km yields a better representation of the concentric pattern of the near-surface fields related to the typhoon center. In general, the use of different flux options results in the enhancement of these surface fields and, thus, the intensity in terms of the minimum central pressure and maximum wind speed, in the order F1 > F2 > F0.

Figure 5Horizontal distributions for (a) wind speed (m s⁻¹) at 10 m height, (b) total rainfall (mm h⁻¹), (c) surface sensible heat flux (W m⁻²), and (d) surface latent heat flux (W m⁻²) at forecast hour 71 h (valid at 22:00 UTC on 7 November) for the experiment F0 with 3 km resolution. All fields are taken from the transformed subset on cylindrical grids; see text for detail. At 71 h, the simulated typhoon centered at the nearest location to the best-track location at 18:00 UTC on 7 November 2013. Range rings are shown at radii of 20, 40, 60, 80, and 100 km. Black dot at the plane center denotes the simulated typhoon center at 71 h, whereas the two dots to the upper left and lower right of this central dot are the corresponding locations of the typhoon center in the next and previous 2 h. Solid lines connecting these dots denote the simulated typhoon track within the 5 h.

Download

Figure 6Horizontal distributions for wind speed at 10 m height (shaded; m s⁻¹), surface latent heat flux (yellow contours; W m⁻²), and total rainfall (cyan contours; mm h⁻¹) at 71 h for all experiments. In all the experiments, the center location of the simulated typhoon at 71 h is at the nearest or the second nearest location to the best-track location at 18:00 UTC on 7 November 2013. The distributions for experiment F0 of 3 km resolution are the same as in Fig. 5 but slightly zoom in on a circle of radius 80 km; range rings are shown at radii of 20, 40, 60, and 80 km. Other symbols are shown as in Fig. 5.

Download

We have examined these surface fields over the period from 12:00 UTC on 7 November to 00:00 UTC on 8 November and found similar features to those in Fig. 6, indicating that this period can be seen as the mature stage of the simulated typhoon. Overall, the horizontal distributions of the wind speed at 10 m height and total rainfall reveal a region of higher wind speed on the right (northern) side of the track and a region of intense rainfall to the left side of the rear side of the track. This result agrees with the study of Shimada et al. (2018) for Typhoon Haiyan. Based on the radar reflectivity measured several hours before the first landfall of Haiyan, they found that the region of the larger maximum wind speed was located on the right side of the track, and the region of stronger retrieved rainfall was located on the downshear-left side. During their analysis period (approximately 2.5 h immediately before Haiyan's first landfall), the 850–200 hPa vertical wind shear was north-northeasterly at about 4 m s⁻¹. Therefore, the downshear-left side is approximately the left-rear side of the track. They found a slight southward vertical tilt (approximately 1 km from 1 to 7 km altitude) of the center of Haiyan, almost the same as the direction of the vertical tilt. We also found similar vertical wind shear of approximately 5 m s⁻¹ at a northeasterly direction, as well as a small vertical tilt (about 1 grid box) of the storm center, during the 12 h period starting at 12:00 UTC on 7 November 2013 (data not shown). Accordingly, our simulations are in reasonable agreement with observations of Haiyan studied by Shimada et al. (2018).

The next step is to examine the azimuthal averages of the radius distribution for surface fields as well as the vertical structure of the simulated typhoon. For the radius distributions shown in Figs. 7, 8, 10, and 11, the track lines over a 3 h window (the two segments centered at that particular hour, as shown in Fig. 5) were used to identify the semicircles on both sides. Because a small deflection of the track tendency is generally found in the 3 h window, the areas of the right (northern) and left (southern) semicircles are not necessarily the same.

Figure 7Radius–time cross sections of (a) 10 m wind speed (m s⁻¹), (b) rainfall (mm h⁻¹), (c) surface sensible heat flux (W m⁻²), and (d) surface latent heat flux (W m⁻²) for the control experiment (3kmF0). The positive and negative distances from the typhoon center represent the azimuthal averages of the right (northern) and left (southern) semicircles, respectively, along the typhoon track. For each hourly output, the corresponding track in 3 h centered at that particular hour is used to identify the semicircles on both sides.

Download

Figure 8The radial variations of surface fields for experiments 1 km (a), 3 km (b), and 6 km (c). The surface fields shown are (from top to bottom) sea level pressure (SLP; hPa), wind speed and tangential and radial winds at 10 m (Wspd, Vt, and Vr; m s⁻¹), surface latent and sensible heat fluxes (LH and SH; W m⁻²), and total rainfall (Rain; mm h⁻¹). Positive values of Vr indicate inflow. Gray, black, and gray-dashed curves represent flux options F0, F1, and F2, respectively. Horizontal gray lines are selected reference lines for each field. All fields are time averaged between 18:00 UTC on 7 November 2013 and 00:00 UTC on 8 November 2013. Positive and negative distances represent the right (northern) and left (southern) semicircles as in Fig. 7. Vertical light-gray and gray shaded bars denote the radial sectors bounded by radii of 50 and 30 km, respectively.

Download

Figure 7 provides Hovmöller diagrams of the 10 m winds, total rainfall rate, and surface enthalpy (latent and sensible heat) fluxes for 3kmF0. Most noteworthy is the asymmetric structure found in the wind and rainfall fields. The stronger winds over the right (northern) side occurred in association with a larger rainfall rate over the left (southern) side of the typhoon center throughout the integration period. The asymmetric structure can also be found in the surface latent heat flux, whereas it is less evident in the surface sensible heat flux. Similar asymmetric patterns are found in all other simulations (data not shown). Figure 8 illustrates the composite of radial distributions of the selected surface fields for the mature stage, namely between 18:00 UTC on 7 November 2013 and 00:00 UTC on 8 November 2013, for the sensitivity experiments F0, F1, and F2. The patterns of F2 experiments are generally somewhat in between those of F0 and F1, with their wind field (enthalpy flux) patterns following those of F1 (F0) more closely. This is expected because F1 and F2 use an identical momentum transfer coefficient and F0 and F2 possess a similar behavior to the enthalpy transfer coefficients in terms of their magnitudes and wind-dependent variations.

We begin with the composite of 3 km resolution (Fig. 8b). During the mature stage, the radial patterns of mean sea level pressure (SLP) and 10 m winds reveal the typical structure of a tropical cyclone. The difference between the radii of maximum 10 m winds (RMWs) of 3kmF0 and 3kmF1 is insignificant. The tangential wind components dominate the 10 m wind fields. In both experiments, the radii of the maximum radial wind component are slightly larger than those of their respective tangential winds. In the context of azimuthal averages, the radial distributions of surface enthalpy fluxes are fundamentally constrained by the 10 m wind speed intensity. The aforementioned asymmetric pattern is again evident in the rainfall distribution and even in the wind fields. The 3kmF1 experiment yields a deeper storm central pressure and stronger wind speeds and enthalpy fluxes than the 3kmF0 experiment does. The 3kmF2 experiment produces similar patterns to the 3kmF0 experiment. The majority of these increments occur around the eyewall (in terms of the maximum winds). Increasing the grid spacing to 6 km yields a composite radial structure highly similar to the 3 km group, thus indicating the relationship among these near-surface fields (Fig. 8c). The major difference between the groups (6 and 3 km) is the size of the inner core with respect to the eyewall radius (measured in RMWs) and the sharpness of peaks in the near-surface fields. This is consistent with the aforementioned horizontal distributions of these surface fields shown in Fig. 6. As grid spacing decreases, RMWs decrease and the peak sharpness in winds and enthalpy fluxes increase. A much more marked decrease in RMWs and increase in peak sharpness are found in the 1 km resolution group than in the other groups (Fig. 8a). The SLPs, wind fields, and enthalpy fluxes in the 1 km group have narrow but sharp peaks, along with larger peak values compared with the other resolution groups. The asymmetric pattern of the wind fields is less evident but still can be seen. The 1kmF1 experiment yields a much stronger storm than the 1kmF0 experiment. The asymmetry of the rainfall pattern, however, decreases by using F1. At this point, whether the reduction of grid spacing or the use of more realistic flux formulas leads to a more symmetric (and thus more intense) tropical cyclone remains unclear.

In brief, both the increased resolution and more reasonable flux formulas enhance storm strength. Storm size (in terms of RMWs) is apparently sensitive to the changes in grid spacing but less sensitive to the choice of flux options. A reduction of storm size with smaller grid spacing has been noted by numerous studies (e.g., Davis et al., 2008; Fierro et al., 2009; Gentry and Lackmann, 2010; Kanada and Wada, 2016). As the strength of a storm increases, the eyewall contracts, indicating a higher efficiency of storm intensification with smaller grid spacing. In our work, the effect of different flux formulas on Typhoon Haiyan simulation is revealed by the changes in storm central pressure and the intensity of maximum winds and enthalpy fluxes around the eyewall. This implies that the underlying mechanisms (i.e., surface flux exchange processes) of various flux options account for storm intensity, whereas the reduction of grid spacing leads to more intense storm structures providing more conducive conditions for enhancing a positive flux effect.

Figure 9 shows the statistics for data points collected within a circle of radius 100 km during a sampling period of 7 h. From the percentages of binned wind speeds, the major pattern of the differences in flux options is revealed by a rightward shift of the curves in the order of F0, F2, and F1. For each resolution group, a reduction of the higher percentages of moderate to high wind speeds and the emergence of higher wind speeds are observed by comparing F2/F1 with F0. This feature is more evident in the distributions of 1 and 3 km resolution groups, in particular for the pattern shift of F1. This is to say that experiments with F1 can produce extremely high wind speeds (exceeds approximately 55 m s⁻¹) as the resolution increases to 3 km as well as 1 km. For the 6 km resolution, less grid points (<5 %) with a wind speed of approximately 55 m s⁻¹ are found even with flux option F1. This can also explain the comparable maximum wind speeds found in the F1 experiments with 3 and 6 km resolution (see Fig. 4c). The comparable storm intensities, in terms of minimum central pressure, found in 3 and 6 km resolutions (see Fig. 4b) are indeed in association with different near-surface wind speed distribution (also see Fig. 6). In 3 km as well as 1 km resolution, more extended and organized areas of extremely high wind speeds are found, whereas for 6 km resolution, broader areas of high wind speeds with a small amount of extremely high wind speed patches are found. The previous analyses indicate that the typhoon intensity is mainly controlled by the model resolution. The near-surface wind speed is an influential factor for the calculation of bulk transfer coefficients. Considering the low frequency of occurrences for extremely high wind speeds, the positive effect of the more reasonable surface option (F1) cannot be enhanced efficiently with a relatively large grid spacing of 6 km.

Figure 9The top panels are the percentage of samples in 2 m s⁻¹ bins for wind speeds. These wind speed samples were collected inside a circle of radius 100 km with respect to the typhoon center at each output hour during a sampling period from 18:00 UTC on 7 November 2013 to 00:00 UTC on 8 November 2013. Over this sampling period of 7 h, only points over ocean are collected. The percentages were the number of data points in each wind speed bin divided by the total count during the sampling period. The upper middle, lower middle, and bottom panels show the binned averages and data ranges for the ratio of C_K∕C_D at 10 m height, latent heat flux, and sensible heat flux, respectively, in each 2 m s⁻¹ wind speed bin. Here the values of C_Q at 10 m height are used as C_K. The binned averages are the mean values for the corresponding variable as collected in bins depending on the collocated wind speed and are displayed as solid or dashed curves. Within each wind speed bin, the respective data range for each variable is bounded by the minimum and maximum values from all binned samples during the period of 7 h. The data ranges are denoted as vertical bars (C_K∕C_D) and symbols (fluxes). Gray curves in the panels of C_K∕C_D are the same as in Fig. 2d, superimposed for reference.

Download

Altogether, the wind speed statistics presented here suggest that more detailed structures over the full wind speed spectrum can only be resolved as the model resolution increases to up to 3 km or higher. This is consistent with the findings of Davis et al. (2008) and Gentry and Lackmann (2010) that a resolution of approximately 3 km or higher is required for significant improvement of intensity simulation of intense hurricanes (Hurricane Katrina in their case). With a lower resolution, 6 km in our case, we speculate that the momentum energy is diluted over the nearby grid points such that limited magnitude (higher but not extreme) of wind speeds would become evenly distributed over a broader area. This idea is specifically inspired by the studies of Bryan et al. (2003) and Miyamoto et al. (2013). We further discuss the energy distribution of these simulations in the next section.

The behaviors of C_K∕C_D show that the patterns almost strictly follow their respective curves as predicted by the calculation of the neutral condition using the pseudo-wind data (see Fig. 2d). The spread of data, as represented by the range, is rather small and is approximately 10 %–12 % of the magnitude of their respective binned averages. However, the resulting latent heat and sensible heat fluxes present a much larger spread of data on the wind speed bins. For lower wind speeds (and thus smaller mean values of fluxes), the range of data is approximately 90 % (150 %) of the magnitude of the binned mean. For high to extremely high wind speeds (and thus large mean values of fluxes), the range of data comes close to or less than 50 % (90 %) of the magnitude of the binned mean. A further inspection of the data reveals that the spread in binned data arises from the discrepancy in stability conditions near the ocean surface. This is expected as seen in Eqs. (A2) and (A3). For moderate to high wind speeds, given a value of latent heat or sensible heat flux, the corresponding wind speed ranges are considerably large. Therefore, the structure of near-surface moist-air enthalpy is also important for typhoon development. Furthermore, F1 can produce the largest latent and sensible heat fluxes between the three flux options, and the magnitude of the interoption differences increases as the resolution increases. This is obviously because of the larger percentage of extremely high wind speeds that F1 can achieve. Although the extremely high wind speeds comprise a small percentage of the eyewall circulation of a mature cyclone, they have a crucial role in enhancing surface fluxes that supply enthalpy to the eyewall.

As the resolution increases, the vertical distribution of tangential winds reveals a significant increase in the vertical extent of stronger wind speed, a reduction of the outward slope of the maximum wind axis, and thus a contraction of the eyewall (Fig. 10). With the change from F0 to F1, tangential wind speeds increase, with an outward and upward expansion of the region of stronger wind speed. No apparent changes are noted in the vertical slope and radius of the eyewall owing to the change in flux options. Both the composite inflow and outflow weaken as the resolution increases but enhance as the flux option changes from F0 to F1. Here, the maximum value of inflow (outflow) is 28.3 (−38.8) m s⁻¹ in 1kmF0 but it is 33.6 (−44.5) m s⁻¹ in 1kmF1. Moreover, the maximum cores of tangential winds descend and the depth of inflow decreases as the resolution increases. Overall, the patterns of F2 are somewhat between those of F0 and F1. In addition, we can see that the upper-level radial inflows are evident in the 1 km resolution group. The upper-level radial inflows are apparently weaker in the 3 km resolution group and are hardly found in the 6 km resolution group. The upper-level radial inflow above the upper outflow layer is related to the development of upper-level warming in the eye; this is addressed in Fig. 12.

Figure 10Radius–height cross sections of tangential wind (contours; m s⁻¹) superposed on radial wind (color shaded; m s⁻¹) for experiments (a) 1kmF0, (b) 3kmF0, (c) 6kmF0, (d) 1kmF1, (e) 3kmF1, (f) 6kmF1, (g) 1kmF2, (h) 3kmF2, and (i) 6kmF2. The tangential winds are contoured every 10 m s⁻¹ starting from 30 m s⁻¹. Positive values of radial wind indicate inflow. The upper-level radial inflows (purple dashed contours; m s⁻¹) are also superimposed and are contoured every 0.5 m s⁻¹ starting from −0.5 m s⁻¹. All fields are time averaged between 18:00 UTC on 7 November 2013 and 00:00 UTC on 8 November 2013. Positive and negative distances represent the right (northern) and left (southern) semicircles as in Fig. 7.

Download

Figure 11The same as Fig. 10 but for vertical velocity (contours; m s⁻¹) and radial wind (color shaded; m s⁻¹). Positive (upward) and negative (downward) velocities are indicated by red and blue contours, respectively, and are plotted every 1 m s⁻¹ starting from wind speeds of 0.5 m s⁻¹. In addition, pink hatched areas denote an upward velocity of 0.3–0.5 m s⁻¹; cyan contours denote a downward velocity of −0.3 m s⁻¹.

Download

Figure 12Time–height cross sections of upper-level radial inflow (contours; m s⁻¹) superposed on the temperature anomaly (color shaded; K) for experiments (a) 1kmF0, (b) 3kmF0, (c) 6kmF0, (d) 1kmF1, (e) 3kmF1, (f) 6kmF1, (g) 1kmF2, (h) 3kmF2, and (i) 6kmF2. The upper-level radial inflows are contoured every 0.5 m s⁻¹ starting from −0.5 m s⁻¹. The temperature anomaly is defined with respect to a composite of domain-averaged temperature profiles at the model initial time of the sensitivity experiments (00:00 UTC on 5 November 2013); the composite of domain-averaged temperature profiles is derived as the mean of eight domain-averaged temperature profiles obtained from all the experiments. The time–height cross section of the temperature anomaly is the area-averaged temperature profile inside a ring of radius 21 km with respect to the typhoon center at each output hour minus the composite of domain-averaged temperature profiles. The upper-level radial inflows are the area-averaged radial winds inside a ring of radius 200 km with respect to the typhoon center at each output hour; only the radial winds above 9 km height are used for calculation.

Download

The contraction of the eyewall with reduction of grid spacing is also revealed by a more upright eyewall updraft (Fig. 11). We start with the major eyewall updrafts and downdrafts as denoted by red and blue contours, respectively, in Fig. 11. The reduction of grid spacing also results in stronger updrafts and downdrafts, where the enhanced eyewall downdrafts may be a cause of the decrease in the depth of inflow. The convections outside of the eyewall (weaker upward velocity denoted by the pink hatched area) are apparent in the 1 km resolution group. For this resolution group, the outside convection in F1 is stronger than that in F0. The outside convection is a result of the azimuthally averaged upward velocity of the outer rainband. In the 3 km resolution group, the outside convections are weaker but still apparently shown. No outside convection can be clearly distinguished from the eyewall in the 6 km resolution group. As the resolution decreases, the model may not be able to resolve intense updraft cores in the eyewall as well as the outer rainband near the eyewall. The larger grid spacing can only resolve a broader area of weaker updraft instead. We further discuss this in the next section. Stronger inflow and outflow as well as broader eyewall updraft are found at the 6 km resolution, suggesting that secondary circulation is stronger at lower resolutions, similar to the finding of Fierro et al. (2009). We further examined the horizontal flow patterns near the altitudes of the upper outflow layer, both on model native grids and cylindrical grids, and confirmed that the upper-level flow patterns in the 6 km resolution group are more radially outward than the rotational component. As the resolution increases, the upper-level flow pattern becomes more rotational, and the radially outward component is suppressed. Because of such discrepancy in the upper-level flow pattern, the azimuthally averaged upper outflows (upper-level radial winds) in the 6 km resolution group are the strongest between the resolution groups.

The time evolution of the core temperature is represented by the temperature anomaly with respect to a reference mean temperature (Fig. 12). Here we term the positive temperature anomaly as warming. The distinct difference found among the resolution groups is the development of the upper-level warming layer within the eye. For the 1 km resolution group, the apparent upper-level warming layer appears above approximately 10 km within a deep warming column in the eye. The upper-level warming layers occur earlier in the experiments F1 and F2 (Fig. 12a, d, g). The warming layer in experiment 1kmF1 is the most intense among all the others. The timings of the upper-layer warming are associated with their respective upper-level radial inflows. The upper-level warming layer has long been recognized by observational studies (e.g., Hawkins and Imbembo, 1976) and numerical simulations (e.g., Zhang and Chen, 2012; Chen and Zhang, 2013; Wang and Wang, 2014). In recent years, the formation of the upper-level warming layer has been related to the rapid intensification of tropical cyclones (e.g., Zhang and Chen, 2012; Chen and Zhang, 2013; Wang and Wang, 2014). This upper-level warming layer, under hydrostatic balance and considering the enhanced effect due to the upper-level dryer and thinner conditions of air, can produce a much greater effect on the surface pressure falls than the lower-level warming (e.g., Zhang and Chen, 2012). Regarding the formation as well as the maintenance of the upper-level warming layer, Zhang and Chen (2012) emphasized the contribution from the high-potential-temperature air detrained from the lower stratosphere, whereas Chen and Zhang (2013) introduced the contribution from the air detrained from the convective bursts (anomalous intense updrafts) in the eyewall. Both studies revealed the importance of the upper-level radial inflow in the warm air detrainment onto the eye region. Wang and Wang (2014) reported that, in an environment with decreasing vertical wind shear, air detrained from the lower stratosphere as well as from the convective bursts could result from the development of convective bursts within the inner-core region.

From experiment F0 to F1, the upper-level warming layer enhances significantly (Fig. 12 a, d). Although both F0 and F1 with 1 km resolution can produce extremely intense updraft cores, the experiment F1 provides a larger energy source (through surface enthalpy fluxes) to the simulated typhoon than experiment F0 does. Therefore, the significantly enhanced upper-level warming layer found in the experiment 1kmF1 could result from the warm detrainment from the intense updrafts (may be seen as convective burst) in the eyewall. The pattern of the experiment F2 is again somewhat between those of F0 and F1. A similar development of the upper-level warming layer can also be found in the 3 km resolution group, but apparently weaker than that in the 1 km group. However, the development of upper-level radial inflow appears to precede the main warming for nearly 24–36 h. This may be partly attributed to the weaker contraction of the eyewall in the 3 km resolution group. The weaker upper-level radial inflow is perhaps another reason for this inconsistency in the time evolution. The 6 km resolution group does not present comparable development of this upper-level warming layer as in the 3 km group, again indicating that a resolution of 6 km is inadequate to resolve the details of structure and development of a strong typhoon like Haiyan (2013).

When the resolution increases, a more intense convection core can be resolved, leading to a rapid pressure decrease near the typhoon and a pressure gradient increase in the area. This change in the pressure field is presumably due to the enhanced diabatic heating near the typhoon eyewall. The larger pressure gradient force would reduce the RMW according to the gradient wind balance equation. The reduced RMW would then accompany an area of diabatic heating closer to the typhoon center. This contraction of the eyewall keeps decreasing the central pressure of the typhoon, yielding positive feedback for typhoon development. The enhancement of these processes is mainly related to the model grid spacings. Near-surface wind speed affects the calculation of bulk transfer coefficients. The positive effect of the more reasonable surface option (F1) cannot be enhanced efficiently unless an extremely large wind speed is generated through the reduction of grid spacing.

Several studies have reported that the reduction of grid spacing yielded a deeper, stronger, and more upright and contracted eyewall (e.g., Fierro et al., 2009; Gentry and Lackmann, 2010; Kanada and Wada, 2016). We also demonstrated that the size of updraft cores in the eyewall shrank and the region of downdraft increases as the resolution is increased. Furthermore, both updraft and downdraft become more intense with the reduction of grid spacing. These features have also been recognized in many studies (e.g., Fierro et al., 2009; Gentry and Lackmann, 2010; Gopalakrishnan et al., 2011; Kanada et al., 2012). In addition, large vertical transport of heat and moisture can result in stronger storms given that the model resolution is sufficiently high (e.g., Gopalakrishnan et al., 2011; Kanada et al., 2012).

We used contoured frequency by altitude diagrams (CFADs; Yuter and Houze, 1995) to examine the variability of convective cores over the region shown in Fig. 3 (the blue box). The CFADs is a statistical method for summarizing the vertical distributions of meteorological fields. A CFAD is constructed by collecting frequency distributions of a particular variable at evenly spaced altitudes within an area, compiling them into a two-dimensional (data bin and altitude) dataset, and portraying the data on a single contour plot. The ordinate in the plot represents the altitude variation, and the abscissa represents the frequency bin. For each altitude on a CFAD, the frequencies should add up to 100 %. For a given CFAD, each point depicts the frequency of occurrence of the data in that bin at a specific altitude. Accordingly, the CFADs ignore horizontal variability and provide a bulk statistical measure for comparing the vertical structure of evolving fields of cumulonimbus clouds or any convective systems. Here, we take Fig. 13a as an example. On this CFAD, there are higher percentages (e.g., bounded by 20 % of data per decibel relative to Z per kilometer, or % dBZ⁻¹ km⁻¹) found in the higher reflectivity bins at lower altitudes and found in the lower reflectivity bins at higher altitudes. The former can indicate convective cells or precipitations, whereas the latter can indicate snow or stratiform precipitation. Therefore, considering a set of CFADs for an evolving cumulonimbus clouds, say from the initiation to the mature stage, the maximum percentage at each altitude would change from vertically oriented to negatively tilted toward lower reflectivity values. More detailed interpretations of the use of CFADs can be found in the work by Yuter and Houze (1995).

Figure 13The CFADs of simulated reflectivity for experiments (a) 1kmF0, (b) 3kmF0, (c) 6kmF0, (d) 1kmF1, (e) 3kmF1, (f) 6kmF1, (g) 1kmF2, (h) 3kmF2, and (i) 6kmF2. The sampling time for model output is taken as the simulated typhoon centered at its nearest location to the best-track location at 18:00 UTC on 7 November 2013 and for the following 2 h. The simulated reflectivity was taken from all grid points with updraft (positive vertical velocity) within the calculation domain for CFADs (see blue box in Fig. 3 for the corresponding area) in the model native grids. For reflectivity CFADs, the bin size is 5 dBZ and the plot is contoured at 0.01 %, 0.1 %, 1 %, 5 %, 10 %, 20 %, 30 %, 40 %, 50 %, and 60 % of data per decibel relative to Z per kilometer (% dBZ⁻¹ km⁻¹), with the 20 % dBZ⁻¹ km⁻¹ contour highlighted in black and the lowest two contours displayed in magenta.

Download

CFADs of simulated reflectivity were conducted separately for updrafts and downdrafts (Figs. 13 and 14). In general, the CFADs of reflectivity derived from both updrafts and downdrafts present similar patterns. In these plots, a region bounded by a specific value, say 20 % dBZ⁻¹ km⁻¹, signifies that 80 % of the sampled data points are within this region. Taking the region bounded by the value of 20 % dBZ⁻¹ km⁻¹ as a reference, these plots reveal that approximately 80 % of the sampled data points in all the simulations present an increase in reflectivity with height above approximately 4 km; at heights lower than 4 km, reflectivity demonstrates little change. Similar patterns can be found in studies of CFADs of reflectivity for strong TCs in numerical simulations (Fierro et al., 2009) and for Florida cumulonimbi revealed by observational radar data (Yuter and Houze, 1995). Most noteworthy about our CFADs of reflectivity for both updrafts and downdrafts is the shrinking region bounded by a contour of 20 % dBZ⁻¹ km⁻¹ as the resolution increases. This indicates that at higher resolutions the range of reflectivity magnitudes presents a relatively broader distribution, which implies an increased occurrence of higher reflectivity. A descending of the core region of higher frequencies (as bounded by a contour of 20 % dBZ⁻¹ km⁻¹) is found in the CFADs of the downdraft compared with those of the updraft. According to the findings of Yuter and Houze (1995), this implies that the larger reflectivity at lower levels may be associated with the downdraft and thus can be interpreted as heavy precipitation, whereas higher reflectivity at approximately 4–5 km is mostly associated with updraft and, thus, may be seen as a convective cell. However, interpretations are different for reflectivity associated with updraft and downdraft here.

Figure 14The same as Fig. 13 but for CFADs of simulated reflectivity with downdraft.

Download

We further conducted CFADs of vertical motions to address this issue (Fig. 15). Most vertical velocities, carried by approximately 95 % of the data points as indicated by the region within the 1 % of data per meter per second per kilometer (% m⁻¹ s km⁻¹) contour, are apparently less than 5 m s⁻¹ in all simulations. A similar pattern of CFADs of vertical velocity, revealing that most common values of vertical velocity are substantially low, can be found in numerical simulations for TCs (Fierro et al., 2009; Gentry and Lackmann, 2010; Wang and Wang 2014) and observational data for cumulonimbi (Yuter and Houze, 1995). Furthermore, the region of upward motion is larger than that of downward motion in general. A similar asymmetry of the distributions of upward and downward motion was evident in simulations in previous studies (Fierro et al., 2009; Gentry and Lackmann, 2010; Wang and Wang, 2014). In Fig. 15, the most significant differences in the CFADs of vertical velocity are revealed by the lower frequencies, namely the values of 1 % m⁻¹ s km⁻¹ and 0.1 % m⁻¹ s km⁻¹. As the resolution increases here, both upward and downward motion present an expansion of the distribution of the frequency of 1 % m⁻¹ s km⁻¹. However, the extents of the changes in upward and downward motions with resolution are different. In all simulations, the changes in the downward motion with an increased resolution are more significant than those in the upward motion. Therefore, the region (occurrence) of the downward motion increases as the resolution increases. The results also indicate the presence of less downdrafts in lower-resolution simulations (i.e., 3 and 6 km), and the degradation of available data points can result in large biases in the CFADs of reflectivity in association with the downdraft.

Figure 15The same as Fig. 14 but for CFADs of simulated vertical velocity. The velocity bin size is 1 m s⁻¹ and the plot is contoured at 0.01 %, 0.1 %, 1 %, 5 %, 10 %, 20 %, 30 %, 40 %, 50 %, and 60 % of data per meter per second per kilometer (% m⁻¹ s km⁻¹), with the 20 % m⁻¹ s km⁻¹ contour highlighted in black and the lowest two contours displayed in magenta.

Download

Finally, the 1 km simulations produce stronger updrafts and downdrafts than the 3 and 6 km simulations. Therefore, the higher reflectivity (indicated by the 0.1 % dBZ⁻¹ km⁻¹ contour) in 1 km simulations may be associated with much more intense updrafts compared with those in the 6 km simulations. This is consistent with the discrepancy between the secondary circulations of different resolutions in the previous section, where 6 km simulations presented relatively wider and weaker updrafts among the different resolutions. As previously mentioned, the CFADs of reflectivity in the 6 km simulations presented a relatively concentrated range of reflectivity magnitudes compared with the 1 km simulations, implying that the convective cells associated with the eyewalls may be broader in spatial scale. This is confirmed by an examination of horizontal distributions of the reflectivity of the 6 km resolution of simulations (data not shown). This can be interpreted as follows. As the resolution is insufficient to resolve the typical cumulonimbus size (∼10 km), a number of convective cells (e.g., cumulonimbus clouds) would be interpreted as a broader and perhaps stronger convective plume as discussed by Bryan et al. (2003) and Miyamoto et al. (2013). Bryan et al. (2003) attributed this broader convective plume to the diffusion due to subgrid-scale processes. In short, the intense reflectivity (convective) cores of the 1 km simulations are driven by extreme updrafts occurring within a rather small fraction of the CFAD's computational domain, whereas the 6 km simulations are driven by relatively weaker updrafts spread over a much broader spatial distribution.

In general, this resolution dependence of the spatial scale of updrafts may be attributed to the model effective resolution, namely the minimum scale resolved by the discretized model, as discussed by Skamarock (2004). By locating where the simulated energy slope drops below the expected slope, Skamarock (2004) found that the effective resolution of the WRF model at a grid spacing of 4–22 km is generally nearly 7-fold that of the model grid spacing (DX) and thus indicated that only information at scales greater than approximately 7 DX in energy spectra represents a physical solution for the WRF model. To examine the effective model resolution in our experiments, we adopted the algorithm of spectral computation introduced by Denis et al. (2002) for the calculation of kinetic energy (KE) spectra. Denis et al. (2002) introduced a two-dimensional discrete cosine transform to convert the limited domain of grid point atmospheric fields into spectral space. In this study, the spectral decompositions of velocities (u, v, and w) were first computed for each hourly model output and for each height level, and the KE spectra were then calculated from these velocity spectra. Finally, the KE spectra were averaged over selected periods and layers for Fig. 16. Because the results of different flux options are similar, only the experiments with F0 and F1 were plotted in Fig. 16. In all simulations, for the KE spectra derived from the horizontal velocity components, we find a model effective resolution of approximately 7 DX where the values of the KE spectra drop significantly, as indicated by Skamarock (2004). Although the vertical velocity spectra exhibit a flatter slope than the horizontal KE spectra, they also present apparent energy drops at wavelengths of approximately 6–7 DX. Our vertical velocity spectra show a pattern (i.e., the slopes) that is somewhat different from that of Bryan et al. (2003) with respect to the absence of peak values around the respective model effective resolutions. Bryan et al. (2003) conducted simulations with horizontal resolutions much higher than 1 km: 500, 250, and 125 m. With their extremely high resolutions, more intense convections could be resolved. They indicated that a model resolution of 1 km remains insufficient for explicitly resolving convective clouds. Here, we speculate that insufficient intense convection cores in our simulations may result in the absence of peak values around the respective model effective resolutions.

Figure 16Mean KE spectra for horizontal wind (a, b) and vertical velocity (c, d) for all simulations as computed over the boundary layer (a, c; between 1 and 2 km) and free troposphere (b, d; between about 3 and 8 km) for the mature stage. The mature stage is indicated by the average between 18:00 UTC on 7 November 2013 and 00:00 UTC on 8 November 2013. The gray lines correspond to the expected slopes for the large scale (−3 slope) and mesoscale to smaller scales ($-\mathrm{5}/\mathrm{3}$ slope) of the atmospheric kinetic energy spectrum. The wavelength and wavenumber are given on the lower and upper x axes, respectively. The purple triangles in each panel, from left to right, denote the locations of 7 DX for the simulations of resolutions 6, 3, and 1 km, respectively. The domain of the spectral decomposition is shown in Fig. 3.

Download

Finally, the slope of the physical portion of the spectra (for the horizontal and vertical components, respectively) remains essentially unchanged as the model resolution is varied. In general, the higher the resolution, the further the downscale extent is, the smaller the scales represented are, and the smaller the model effective resolutions reached are. The aforementioned descriptions remain true for all flux options and averaging layers in our work. The key point is that the effective resolution is determined by grid spacing, not the flux options. We noted some subtle differences between the F0 and F1 experiments at smaller wavelengths at nearly 2 DX in the KE spectra but in only a few hours of the entire simulation period. In short, although the use of more reasonable flux formulas can increase simulated storm intensity to some extent, the positive effect of the flux formulas cannot be efficiently enhanced unless the grid spacing is appropriately reduced to yield intense and contracted eyewall structure.

The National Centers for Environmental Prediction Global Forecast System (GFS) operational global analysis dataset can be downloaded from the website https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs/ (last access: 20 July 2019; Environmental Modeling Center 2003, 2016). The RTG-SST dataset is available for download at ftp://polar.ncep.noaa.gov/pub/history/sst/ (last access: 20 July 2019; Gemmill et al., 2007). The WMO subset of the IBTrACS (IBTrACS-WMO, v03r09) can be accessed at ftp://eclipse.ncdc.noaa.gov/pub/ibtracs/v03r09/wmo/ (last access: 20 July 2019; Knapp et al., 2010). Appendix B provides the information of the observational estimates of surface flux.

The supplement related to this article is available online at: https://doi.org/10.5194/nhess-19-1509-2019-supplement.

MTK designed the experiments and analyzed and prepared the manuscript with contributions from all co-authors. WMC, YFS, SCL, and EY performed the simulations and processed and analyzed the model output. TRW, YLT, and CYL joined the discussions of the study. All authors contributed to the preparation of this paper.

The authors declare that they have no conflict of interest.

This research has been supported by the Ministry of Science and Technology, Taiwan (grant nos. 105-2111-M-001-003 and 106-2111-M-001-009).

This paper was edited by Ricardo Trigo and reviewed by two anonymous referees.