Numerical Study of Boundary Layer Structure and Rainfall after Landfall of Typhoon Fitow (2013): Sensitivity to Planetary Boundary Layer Parameterization

The planetary boundary layer (PBL) determines the momentum, heat and water vapor exchanges, and vertical mixing between the underlying surface and the atmosphere above (Chen and Ding, 1979; Stull, 1988). Previous studies have shown that surface fluxes and vertical mixing in the boundary layer play a key role in determining changes in the intensity and structure of tropical cyclones (TCs) (Malkus, 1958; Malkus and Riehl, 1960; Ooyama, 1969; Rosenthal, 1971; Emanuel, 1986, Emanuel, 1995, Emanuel, 1997; Rotunno and Emanuel, 1987; Wang et al., 2001; Davis and Bosart, 2002; Deng et al., 2005). Some significant field experiments, such as CBLAST (Coupled Boundary Layer Air-Sea Transfer; Black et al., 2007), have demonstrated that PBL processes play an important role in air-sea interactions during the evolution of TC intensity and structure. Meanwhile, the CBLAST experiment has also led to many important new findings in the rarely sampled hurricane boundary layer. Additional studies (Hill and Lackmann, 2009; Foster, 2009; Zhu et al., 2014) have further confirmed the importance of eddy diffusivity in simulations of TC intensity and structure. For instance, (Zhu et al., 2014) demonstrated that the vertical turbulent mixing scheme played a significant role in asymmetric structures and eyewall mesovortices in TCs. The eddy exchange coefficients are the key factors regulating the characteristics of eyewall disturbances. On the basis of new observations, vertical eddy diffusivity parameterization was improved and tested in both idealized and real hurricane simulations by a hurricane prediction model (Gopalakrishnan et al., 2013; Zhang et al., 2015). These studies illustrated substantial improvements in TC size and surface inflow angle. Recently, the effects of vertical eddy diffusivity parameterization, as well as surface flux parameterization, on the evolution of landfalling hurricanes have been investigated (Ming and Zhang, 2016; Zhang and Pu, 2017; Zhang et al., 2017). The results showed that the modification in the K_m parameter (the vertical diffusion coefficients of momentum) and strong vertical mixing lead to improved simulations of hurricane track, intensity and quantitative precipitation.

Studies have also shown that TC intensity, structure, and precipitation are very sensitive to the choice of PBL parameterization scheme. For example, (Braun and Tao, 2000) examined the sensitivity of high-resolution simulations of Hurricane Bob (1991) to PBL parameterization schemes [including the bulk aerodynamic scheme, Blackadar scheme, Medium-Range Forecast (MRF) scheme, and Burk-Thompson scheme] using the fifth-generation Pennsylvania State University-National Center for Atmospheric Research Mesoscale Model. These authors found that the Burk-Thompson and bulk aerodynamic schemes produced the strongest storms, while the MRF scheme produced the weakest storm. Moreover, the simulated distribution of precipitation also varied substantially among different PBL schemes in a study by (Braun and Tao, 2000). (Li and Pu, 2008) performed simulations using different PBL schemes and showed that the differences in the mean sea level pressure between the Yonsei University (YSU) and Mellor-Yamada-Janjić (MYJ) schemes were up to 19 hPa during the early rapid intensification of Hurricane Emily (2005). (Smith and Thomsen, 2010) compared several PBL schemes (including local, non-local, and higher-order closure schemes) in idealized f-plane simulations and found significant variations in storm structure, intensity, and intensification rate, even when the same surface-flux parameterization scheme was used in all simulations. Several PBL parameterizations were reviewed and assessed by (Kepert, 2012) using a diagnostic TC model. The results indicated that the Louis scheme and a higher-order closure scheme were suitable and recommended for TC simulation. (Wang, 2013) studied the sensitivity of the simulated Typhoon Muifa (2011) to PBL schemes and showed that both TC intensity and track were sensitive to the PBL scheme, with the YSU scheme producing a more realistic intensity and track for a weak TC than the MRF scheme. (Liu et al., 2017) revealed that the stronger surface fluxes and vertical mixing in the MYJ PBL scheme than in the YSU PBL scheme lead to enhanced air-sea interaction, which helps generate more realistic simulations of the rapid intensification process of Hurricane Katrina (2005) offshore before landfall.

The sensitivity of the simulated TC intensity and structure to the choice of PBL scheme stems primarily from different assumptions adopted in different PBL schemes in considering the transfer of mass, moisture, and energy, which may lead to differences not only in the boundary layer but also in the whole troposphere (Zhang and Zheng, 2004; Hu et al., 2010; Zhang et al., 2011; Wang, 2012; Xu et al., 2013). Although previous studies have illustrated the sensitivity of TC simulation and prediction to PBL parameterization schemes, most of those studies focused on idealized TC simulations or offshore TCs. So far, little attention has been given to the potential sensitivities of the boundary layer structure and precipitation in landfalling TCs (which commonly produce tremendous damage in coastal areas) to PBL parameterization in numerical models.

The objectives of this study are two-fold: to analyze the simulated boundary layer structure and precipitation in the landfalling Typhoon Fitow (2013) and to provide a preliminary study to examine the sensitivities of the simulations to different PBL parameterizations using the WRF model. The rest of the paper is organized as follows. Section 2 gives an overview of Typhoon Fitow (2013). The numerical experiments and data used in this study are described in section 3. The results from the control and other sensitivity experiments are analyzed in Section 4. The main conclusions are summarized in the last section.

As shown in Fig. 1a, Typhoon Fitow (2013) was detected as a tropical storm over the sea east of the Philippines in the western North Pacific on 30 September 2013. Typhoon Fitow (2013) then moved northward and intensified to a severe typhoon at 0900 UTC 4 October 2013. Fitow turned northwestward as a severe typhoon and made landfall over the southeast coast of China at 1715 UTC 6 October 2013, with a central sea level pressure of 955 hPa and a maximum sustained near-surface wind speed of 42 m s^-1. After landfall, Fitow moved westward and weakened quickly to a tropical depression at 0100 UTC 7 October 2013, and finally dissipated over land in Fujian Province (denoted as FJ in Fig. 1a). Fitow was the strongest typhoon that made landfall in China in October from 1949 to 2013 and caused severe damage to Zhejiang Province (denoted as ZJ in Fig. 1a) due to its strong winds and heavy rainfall. The maximum wind gust speed of 76.1 m s^-1 was observed at Shipingshan station (on the southeast coast of Zhejiang Province), which set a new record for TCs over China. Fitow produced the strongest large-area daily precipitation in Zhejiang Province (100 000 km²), reaching 137.5 mm on 7 October 2013. Consequently, Fitow resulted in severe urban water logging, flooding, and debris flow disasters, and a direct economic loss of 27.6 billion Yuan RMB (approximately 4.5 billion US dollars) in Zhejiang Province.

Figure 1.  (a) Track of Typhoon Fitow (2013). (b) Radar reflectivity mosaic captured by the Zhejiang Doppler radar network at 1600 UTC 6 October 2013. (c, d) Synoptic weather patterns at (c) 1800 UTC 6 October 2013 and (d) 0600 UTC 7 October 2013, from NCEP reanalysis data, showing geopotential height (units: gpm; solid contours with intervals of 40 gpm), wind (vectors; units: m s^-1), and relative vorticity (shading; units: 10^-5 s^-1) at 850 hPa. The black solid arrow in (b) indicates the vertical wind shear.

Figure 1b shows the composite radar reflectivity, defined as the maximum reflectivity in each grid volume, captured by the Zhejiang Doppler radar network at 1600 UTC 7 October 2013, when Fitow was about to make landfall accompanied by heavy precipitation. The radar reflectivity displays a significant asymmetric distribution with strong echoes (≥ 30 dBZ) to the north of the storm center, covering almost the entirety of Zhejiang Province. Additionally, the field of vorticity illustrates a clear asymmetric pattern, with the strong cyclonic vorticity mainly in the northern part of Fitow (Fig. 1c). The asymmetry was most likely related to the strong southwesterly vertical wind shear (indicated by the black solid arrow in Fig. 1b) and the westward translation of the storm, in line with the results of (Yu et al., 2015). Strong convection tended to appear on the downshear side and to the right of the typhoon track. Rainfall during the first 12 h of landfall (from 1200 UTC 6 October 2013 to 0000 UTC 7 October 2013) was the strongest and was concentrated in the coastal region of Zhejiang Province. As Fitow weakened quickly over land, rainfall during the second 12 h (from 0000 to 1200 UTC 7 October 2013) was primarily located in northern Zhejiang Province, which resulted mainly from the interaction between an inverted V-shaped trough of the remnant vortex and the easterly jet sustained by TC Danas (2013) to the east of Fitow (Fig. 1d).

The WRF model, version 3.4.1 (Skamarock et al., 2008), was used to conduct a 6-h dynamical initialization (from 0000 to 0600 UTC 6 October) and 30-h simulations of Typhoon Fitow (2013) (from 0600 UTC 6 October 2013 to 1200 UTC 7 October 2013) with triple-nested domains (see Fig. 2a). The outermost domain had 271× 211 grid points with an 18-km horizontal grid spacing. The sizes of the intermediate domain (D02) and innermost domain (D03) were 421× 391 grid points with 6-km horizontal grid spacing and 571× 601 grid points with 2-km horizontal grid spacing, respectively. There were 36 full terrain-following σ levels in the vertical direction for all domains, with 13 levels below 700 hPa and a model top at 50 hPa. The time steps were 90, 30 and 10 s for D01, D02 and D03, respectively. The dynamic initialization scheme developed by (Cha and Wang, 2013) was used for the outer two domains to improve the initial condition of the TC vortex. Specifically, the dynamic initialization span up the axisymmetric inner core of the typhoon vortex by cycle runs in a 6-h window from 0000 to 0600 UTC 6 October 2013 [refer to (Cha and Wang, 2013) for details]. Namely, the outer two domains (D01 and D02) were both started at 0000 UTC 6 October 2013, while the innermost domain (D03) was activated at 1200 UTC 6 October 2013, and integrated over 24 h from 1200 UTC 6 October to 1200 UTC 7 October 2013. Following (Wang, 2013), large-scale spectral nudging was applied in the outermost domain (D01) during both the model spin-up and the subsequent model integration in all experiments. This nudging was performed to eliminate possible bias in the simulated large-scale environmental flow, because our main focus in this study is on the sensitivity of the boundary layer structure and rainfall associated with Typhoon Fitow (2013) to the model's PBL parameterization. The model physics included the WRF single-moment 6-class microphysics scheme (Hong et al., 2004), the Dudhia shortwave radiation scheme (Dudhia, 1989), the Rapid Radiative Transfer Model for longwave radiation calculations (Mlawer et al., 1997), the revised Monin-Obukhov surface layer scheme (Jiménez et al., 2012), and the Noah land surface model (Chen and Dudhia, 2001). The Betts-Miller-Janjić cumulus parameterization scheme (Janjić, 1994, 2000) was applied only in D01.

Figure 2.  (a) Configuration of the model domain. (b, c) Fitow's (b) tracks and (c) maximum wind speeds (units: m s^-1) from 0000 UTC 6 October 2013 to 0000 UTC 7 October 2013, from observation and simulations. The outermost domain, D01, in (a) corresponds to the 18-km horizontal resolution, while domain D02 (indicated by the white frame) has a 6-km resolution, and domain D03 (indicated by the black frame) has a 2-km resolution. The dashed lines in (b) indicate the 12-h simulated tracks from 0000 to 1200 UTC 7 October 2013.

There are two categories for the PBL parameterization in the WRF model. The first category is the so-called first-order closure, which includes the MRF scheme (Hong and Pan, 1996) and the YSU scheme (Hong et al., 2006). The YSU PBL scheme is a revised version of the MRF scheme with the inclusion of an explicit treatment of entrainment processes at the top of the PBL. The second category is higher-order closure, which includes the Bougeault-Lacarrere (BouLac, Bougeault and Lacarrere, 1989), MYJ (Janjić, 2001), and Mellor-Yamada-Nakanishi-Niino Level 2.5 (MYNN2; Mellor and Yamada, 1982, Nakanishi and Niino, 2004) schemes. The BouLac scheme is a 1.5-order closure scheme with the characteristic length for turbulent eddies recomputed by limiting the length to be within a certain range. The MYJ scheme represents a nonsingular implementation of the Mellor-Yamada level 2.5 turbulence closure model through the full range of atmospheric turbulent regimes. MYNN2 is also a turbulent kinetic energy (TKE) closure scheme, based on the Mellor-Yamada model (Nakanishi and Niino, 2004), with both level 2.5 and level 3.0 options (the 2.5-level scheme is used in our study). In addition, the Quasi-Normal Scale Elimination (QNSE) scheme (Sukoriansky et al., 2005; Sukoriansky and Galperin, 2008) uses higher-order closure in unstable stratification and a K-ε model that employs the vertical eddy viscosities and eddy diffusivities developed from turbulence spectrum closure in stable stratification.

In this study, two groups of experiments were performed for Typhoon Fitow (2013) (Table 1). The first group (referred to hereafter as group 1) included the YSU, MYNN2 and BouLac schemes, and featured the same surface layer. The second group (referred to hereafter as group 2) included the MYJ and QNSE schemes, and featured different surface layers. Considering the good general performance of the simulation with the YSU scheme with regard to the track, intensity (see section 4.1, Table 2, Fig. 2b, and Fig. 2c), and structure of Fitow (see section 4.2 and Fig. 3), as well as the related precipitation (see section 4.3), for convenience, we treat this simulation as the control run in the following discussion. The initial and lateral boundary conditions of the model in all experiments were taken from the Global Forecasting System FNL (final) analysis data of the National Oceanic and Atmospheric Administration (NOAA) National Operational Model Archive and Distribution System, which has a horizontal resolution of 0.5^° latitude by 0.5^° longitude on 26 pressure levels at 6-h intervals.

Figure 3.  The composite radar reflectivity (shaded; units: dBZ) from (a) observation and (b-f) all simulations [(b) YSU, (c) MYNN2, (d) BouLac, (e) MYJ and (f) QNSE] at (a1-f1) 1200 UTC, (a2-f2) 1500 UTC, (a3-f3) 1800 UTC, and (a4-f4) 2100 UTC 6 October 2013.

The best track data for Typhoon Fitow (2013) were obtained from the Shanghai Typhoon Institute of the China Meteorological Administration (Ying et al., 2014). The hourly precipitation observations were obtained from approximately 2000 stations of the Zhejiang Automatic Weather Station (AWS) network, and interpolated to the model grids using the inverse distance squared weighting method. According to our own tests and those of other studies (e.g., Ikeda et al., 2010), the result was insensitive to the interpolation method for data with a high resolution (approximately 8 km in this study). The analyzed precipitation data had an hourly temporal resolution and covered the area (26^°-32^°N, 117^°-123^°E).

Considering the comparability between the best-track data and the simulation, we obtained the track and intensity from simulation D02. Figure 2b shows the 24-h best track from 0000 UTC 6 October 2013 to 0000 UTC 7 October 2013, and the 36-h simulated tracks from all five experiments with an additional 12 h from 0000 to 1200 UTC 7 October (denoted with dashed curves) for Typhoon Fitow (2013). The additional 12-h tracks are used in the track sensitivity analysis (see Table 2), as well as for precipitation in section 4.3. All experiments reproduced the storm track reasonably well, especially in group 1. Overall, for the 18-h mean error from 0600 UTC 6 October 2013 to 0000 UTC 7 October 2013, illustrated in the second column of Table 2, the BouLac and YSU schemes performed quite well among the five PBL schemes, with a position error as low as approximately 19 km. The maximum track error was only 42 km in the simulation with the QNSE scheme. This discrepancy became larger after landfall, with the most notable northern track coming from the YSU scheme and the most notable southern track from the QNSE scheme. If we use the YSU simulation (control run) as the baseline, the mean distance after landfall between the QNSE simulation and the baseline simulation from 1800 UTC 6 October 2013 to 1200 UTC 7 October 2013, reached 87 km (see the fourth column of Table 2), which was nearly three times that between the BouLac simulation and the baseline simulation. Specifically, at 1200 UTC 7 October 2013, the difference in position distance between the simulations with the QNSE scheme and the YSU scheme was as large as 115 km.

Because TC motion is mainly controlled by the large-scale steering flow, we calculated the steering flow and its zonal and meridional components from 0600 UTC 6 October 2013 to 1200 UTC 7 October 2013 for all experiments, as listed in Table 3. Note that the steering flow was defined as the mean vector wind averaged within a radius of 500 km from the storm center between 300 and 850 hPa (Holland, 1984). As expected, the southward steering flow of 0.76 m s^-1 was in the simulation with the QNSE scheme, while the northward steering flow of 0.17 m s^-1 was in the simulation with the YSU scheme. This finding explains well the difference in the simulated storm tracks and suggests that the PBL scheme may affect the steering flow and lead to deflection in the simulated track. The track difference also partly explained the difference in the simulated precipitation distribution, as will be discussed in section 4.3.

Figure 2c shows the evolution of the observed and simulated storm intensity. All simulations reasonably captured the intensity change, including an initial 6-h spin-up and the rapid weakening after landfall, although the weakening was slower than that observed. Among all tests, the minimum mean absolute error (MAE) in the simulated storm intensity was only 2.8 m s^-1 with the YSU scheme in group 1, while the maximum error was 4.1 m s^-1 in the group 2 simulation with the QNSE scheme (listed in Table 2). The QNSE scheme produced the weakest storm, with a peak intensity of 40 m s^-1, and the BouLac scheme produced the strongest storm, with a peak intensity of 45 m s^-1. Moreover, the differences among all simulations tended to decrease after Fitow made landfall, partly due to the rapid weakening of the storm. The differences in the simulated intensity among the five PBL schemes were generally less than 6 m s^-1, which indicates that the PBL scheme affected the simulated storm intensity for Fitow, but the difference was not as large as that reported in some earlier studies, as mentioned in section 1. The short maintenance period of Fitow over land, as is typically observed for most landfalling TCs, was probably a reason leading to this difference.

Considering more detailed features in the simulations with high resolution, we took the results from D03 to carry out the analyses of Fitow's structure and precipitation. First, we verified the simulated reflectivity to ensure the reliability of the fundamental structure of Fitow for further analysis. Figure 3 compares the composite reflectivity at every 3-h interval from all simulations with that composited from six Doppler radars in Zhejiang Province. All simulations reproduced the structural evolution of the Fitow landfall process. Prior to landfall at 1200 UTC 6 October, Fitow had a complete eyewall structure (Figs. 3a1-f1). Note that the broken eyewall to the south of the storm center in the observation was due to the missing data (beyond the radar's maximum detection range). Strong spiral rainbands with echoes greater than 30 dBZ developed in the coastal area of Zhejiang Province. However, the eye geometry of Fitow varied with different PBL schemes. Group 1 roughly produced an elliptical eye with a long axis in the east-west direction (Figs. 3b1-d1, Figs. 3b2-d2); similarly, QNSE simulated an elliptical eye with a long axis in the east-west direction (Figs. 3f1 and f2), while the MYJ scheme produced a quasi-circular eye (Figs. 3e1 and e2). Meanwhile, the sizes of the eye in group 1 were larger than those in group 2. As Fitow approached land and decayed, the eyewall gradually collapsed, while strong spiral rainbands in the northern part of the TC were sustained (Figs. 3a2-f2, 3a3-f3). By 2100 UTC 6 October 2013 (Figs. 3a4-f4), after landfall, the two main rainbands in the southern and northern parts of Zhejiang Province were reasonably simulated overall. In comparison, the rainband intensity and location differed considerably with various PBL schemes.

Figure 4 shows the vertical cross sections of the azimuthal mean tangential (shading) and radial (contours) winds averaged during the 6-h period from 1200 to 1800 UTC 6 October 2013, (Figs. 4a1-e1) and from 1800 UTC 6 October 2013 to 0000 UTC 7 October 2013 (Figs. 4a2-e2) in all simulations. Overall, all simulations produced the maximum tangential wind at a height of approximately 1.2 km during the first 6 h, approximately 70 km away from the storm center in horizontal distance. Meanwhile, the tangential wind decreased with time. The major differences in the tangential winds among all simulations were found in the maximum and the radial extent of strong winds. The tangential winds were stronger in the simulations with the YSU and BouLac schemes than those in the other three simulations. If the width of strong tangential winds at 1.2 km was defined as the radial extent of tangential wind speeds greater than 24 m s^-1, the strong tangential wind width of 240 km in the simulation with the QNSE scheme was the largest, while that of 210 km in the simulation with the BouLac scheme was the smallest, and that of approximately 225 km in the other three simulations was intermediate. The tangential wind speed in all simulations decreased remarkably in the following 6 h, but the difference in the width among the five simulations remained similar.

Figure 4.  Vertical-radius cross sections of the azimuthal mean tangential (shaded; units: m s^-1) and radial (contours with an interval of 2 m s^-1 from -16 m s^-1 to 16 m s^-1; units: m s^-1) winds averaged across the 6-h periods from (a1-e1) 1200 to 1800 UTC 6 October 2013 and (a2-e2) from 1800 UTC 6 October 2013 to 0000 UTC 7 October 2013, for all simulations: (a) YSU; (b) MYNN2; (c) BouLac; (d) MYJ; and (e) QNSE.

Furthermore, the profile of tangential wind speed (TWS) was an important indicator of TC boundary characteristics (Kepert, 2012). To demonstrate the total TWS sensitivity to the PBL scheme, the 6-h mean TWS vertical profiles averaged within a radius of 300 km from the storm center were compared from 1200 to 1800 UTC 6 October 2013 (Fig. 5a) and from 1800 UTC 6 October to 0000 UTC 7 October 2013 (Fig. 5b). Clearly, the height of the maximum TWS in all simulations was located between 1.0 km and 1.5 km (Figs. 5a and b). The main differences were in the strength at its maximum level (approximately 1.2 km) and that near the surface. The QNSE scheme produced the largest maximum TWS, while the BouLac scheme produced the smallest TWS, and the other three schemes produced approximately similar intermediate TWS values. The largest variation of 2.5 m s^-1 accounted for 11% of the TWS in the simulation with the YSU scheme. The highest vertical shear of tangential wind below a height of 1.5 km appeared in the simulation with the QNSE scheme, and the second highest was that with the MYJ scheme, followed by the MYNN2, YSU and BouLac schemes, in that order.

Figure 5.  Vertical profiles of the 6-h mean area-averaged (a, b) tangential wind (units: m s^-1), (c, d) radial wind (units: m s^-1), (e, f) potential temperature (units: K), and (g, h) water vapor (units: g kg^-1) within a radius of 300 km from the storm center from 1200 to 1800 UTC 6 October 2013 (left-hand panels) and from 1800 UTC 6 October 2013 to 0000 UTC 7 October 2013 (right-hand panels), for all simulations. Note that the vertical scale is 3.0 km in (a-d) and 2.0 km in (e-h).

The radial wind in the vertical-radius distribution showed the maximum inflow near the surface and approximately 100 km away from the storm center (Fig. 4). The inflow gradually decreased with height and became outflow from approximately 1.8 km. The major differences among all the simulations were the strength and depth of the inflow boundary layer. From 1200 UTC to 1800 UTC 6 October 2013 (Figs. 4a1-e1), the inflow strength varied from -15 m s^-1 to -12 m s^-1 in the QNSE, MYJ, YSU, MYNN2, and BouLac schemes, in order, and there was stronger inflow in group 2 than in group 1. During the following 6 h (Figs. 4a2-e2), the inflow strength in all simulations decreased gradually while maintaining the same order. Simulations with the YSU and BouLac schemes produced a deeper inflow layer of up to 2.1 km, while the inflow layer in the other three experiments mainly extended only to approximately 1.7 km in height.

Additionally, the evolution of radial wind profiles was examined (Figs. 5c and d). In the first 6 h of simulation from 1200 to 1800 UTC 6 October 2013 (Fig. 5c), the largest difference in the radial wind between group 1 and group 2 reached 3.0 m s^-1 at 0.3 km, with the strongest inflow of 10.5 m s^-1 in the simulation with the QNSE scheme and the weakest inflow of 7.5 m s^-1 in the simulation with the BouLac scheme. The difference accounted for approximately one-third of the inflow of 9 m s^-1 near the surface in the YSU scheme. The inflow weakened in the following 6 h in all simulations and varied from 5.5 m s^-1 to 7.5 m s^-1 (Fig. 5d). In particular, the inflow layer was approximately 400 m deeper in the simulations with the YSU and BouLac schemes than in the other three schemes.

Meanwhile, investigation of the simulated potential temperature and water vapor mixing ratio revealed that the main difference occurred in the lower boundary layer below 1.2 km (Figs. 5e-h). The YSU and BouLac schemes produced drier and warmer boundary layers, while the QNSE and MYJ schemes produced wetter and colder boundary layers. The maximum differences in the potential temperature and water vapor mixing ratio among the five simulations reached approximately 1.0 K and 0.5 g kg^-1, respectively. Additionally, group 1 simulated a very similar vertical distribution. The largest difference appeared between the simulations with different surface layer schemes (namely, between YSU and QNSE). This difference seems to suggest that in addition to the vertical diffusivity, the algorithm for surface stress and flux calculations is critical to the simulations of the boundary layer structure of TCs as well.

To understand the possible reasons for the differences discussed above, we investigated the spatial distribution and temporal evolution of the vertical mixing and surface fluxes for all simulations. As shown in Fig. 6, the vertical diffusion coefficients of momentum (K_m) and heat (K_h) varied greatly among the five PBL schemes. From 1200 UTC 6 October 2013 to 0000 UTC 7 October 2013, on average, the vertical mixing was markedly stronger in group 1 than in group 2. The maximum vertical diffusion coefficient in the BouLac scheme was more than three times that in the QNSE scheme. Specifically, the height with the maximum vertical diffusion coefficient was higher in the simulation with the YSU scheme than in other simulations; namely, approximately 1.2 km versus 0.4 km in the first 6-h simulation and approximately 1.2 km versus 0.8 km in the later 6-h simulation. Meanwhile, the magnitude of the vertical diffusion coefficient in the YSU scheme showed little change with time, while those in other schemes decreased with time. The different evolution of the YSU scheme compared to the other schemes was related to the different assumptions and algorithms in various schemes (Shin and Hong, 2011; Kepert, 2012). The K_m values in the YSU scheme, as a first-order and non-local scheme, are mainly determined by wind speed and virtual potential temperature structure from the surface up to the top of the PBL (Hong and Pan, 1996; Hong et al., 2006). In contrast, in the other four higher-order and local schemes, the K_m values are basically determined by turbulent energy and thus vertical wind shear at each level, and depend upon conditions at that level (Mellor and Yamada, 1982; Bougeault and Lacarrere, 1989; Janjić, 2001; Nakanishi and Niino, 2004; Sukoriansky et al., 2005; Sukoriansky and Galperin, 2008). These differences contributed to a lower sensitivity to the decay of Fitow after landfall in the YSU experiment (Kepert, 2012). Meantime, the higher altitude of maximum wind speed in the YSU experiment than that of the vertical wind shear apexes in the other four experiments may also have potentially contributed to the greater height of maximum K_m/K_h in the YSU experiment. Consequently, the stronger vertical mixing in group 1 more efficiently transferred moisture and energy upward, which resulted in relatively smaller vertical variations in the boundary layer structure in terms of dynamics and thermodynamics (as displayed in Fig. 5).

Figure 6.  As in Fig. 4 but for the momentum diffusion coefficients K_m (shaded; units: m² s^-1) and heat diffusion coefficient K_h (contours with 10 m² s^-1 intervals from 10 m² s^-1 to 180 m² s^-1; units: m² s^-1).

In general, comparison of the surface fluxes indicated stronger moisture fluxes and weaker sensible heat fluxes (Fig. 7) in QNSE than in group 1. For instance, at 1500 UTC 6 October 2013 (Figs. a1-e1, Figs. f1-j1), the maximum moisture fluxes in the Fitow eyewall for QNSE were larger than 47× 10^-5 kg m^-2 s^-1, while those in group 1 were approximately 39× 10^-5 kg m^-2 s^-1. Additionally, the averaged moisture (sensible heat) fluxes within a radius of 300 km from the TC center in QNSE were stronger (weaker) than those in group 1. These discrepancies——namely, stronger moisture fluxes and weaker sensible heat fluxes in QNSE——may partly explain the somewhat higher water vapor and lower potential temperature near the surface.

Figure 7.  (a-e) Surface moisture fluxes (units: 10^-5 kg m^-2 s^-1) and (f-j) sensible fluxes (units: 10 W m^-2 s^-1) at (a1-j1) 1500 UTC and (a2-j2) 2100 UTC 6 October 2013, from all simulations: (a, f) YSU; (b, g) MYNN2; (c, h) BouLac; (d, i) MYJ; and (e, j) QNSE.

Figure 8 shows the distribution of the observed and simulated 24-h accumulated precipitation from 1200 UTC 6 October 2013 to 1200 UTC 7 October 2013. The simulation with the YSU scheme reproduced the extreme precipitation induced by Fitow in Zhejiang Province reasonably well. The simulated areal mean rainfall was 137.6 mm, which was close to the 137.5 mm in observations. The mean error (ME), root-mean-square error (RMSE), and MAE were 0.16 mm, 67.14 mm, and 45.07 mm (Table 4), respectively. Moreover, the spatial correlation coefficient between the simulation and the observation reached 0.66. The threat score (TS) and equitable threat score (ETS) were as high as 0.89 (0.76) and 0.40 (0.46) for the threshold of 50 (100) mm, respectively. In particular, the simulation with the YSU scheme successfully reproduced the rainfall in northern Zhejiang Province. All four other simulations also captured the main features of the rainstorm in the coastal region and northern Zhejiang Province, but with relatively lower scores. Note that group 1 simulated stronger precipitation in Zhejiang Province than group 2. Furthermore, the simulated precipitation patterns in group 1 with the same surface layer scheme were more similar to each other than those in group 2 with different surface layer schemes. The main precipitation difference among all simulations appeared in the northern area, as indicated by the thick dashed box in Fig. 8, which will be referred to as the study region (SR) in the following discussion.

Figure 8.  Observed and simulated 24-h accumulated precipitation (units: mm) from 1200 UTC 6 October 2013 to 1200 UTC 7 October 2013: (a) observation, and (b-f) simulations with the (b) YSU, (c) MYNN2, (d) BouLac, (e) MYJ, and (f) QNSE PBL schemes. The dashed rectangle (30^°-31^°N, 119.2^°-121^°E) in each panel denotes the region, referred to as the study region (SR), wherein the major precipitation difference occurred among all experiments.

The simulation with the YSU scheme also captured the rainfall evolution in observations reasonably well, as evidenced by the 6-h accumulated rainfall (Fig. 9). The averaged correlation coefficient of the 6-h accumulated rainfall reached up to 0.63 (Fig. 9b), and the averaged TS and ETS for heavy rainfall (more than 13 mm in 6 h) were 0.71 and 0.36, respectively (Figs. 9c and d), suggesting that the simulation with the YSU scheme was skillful for this case. The major difference in rainfall between the YSU simulations and other simulations in the SR occurred in the following 12 h, from 0000 to 1200 UTC 7 October 2013 (Fig. 9). The RMSEs in group 2 were substantially larger (Fig. 9a), and the spatial correlation coefficients, TSs and ETSs were markedly lower than those in group 1 (Figs. 9b-d). For example, the TSs and ETSs were only approximately half of those in the simulation with the YSU scheme (Figs. 9c and d). The results strongly suggest that both the amount and distribution of rainfall were very sensitive to the choice of PBL scheme and surface layer scheme.

Figure 9.  Comparison of the (a) RMSE (units: mm), (b) correlation coefficient (CC), (c) TS, and (d) ETS for the simulated accumulated precipitation for every 6-h period from 1200 UTC 6 October 2013 to 1200 UTC 7 October 2013, among the five experiments. Heavy rainfall is verified with the threshold of 25 mm in [panels (c) and (d)].

Furthermore, in terms of hourly rainfall, the major difference in the SR was found during 0600 to 0900 UTC 7 October 2013 (Fig. 10). The MYJ and YSU schemes produced the weakest and strongest hourly rainfall, with a bias of approximately 5 mm. The heavy rainfall was located in the central and western parts of the SR in the simulation with the YSU scheme and in the northwestern part in the simulation with the QNSE scheme. The simulated rainfall exhibited a scattered distributed in the other three simulations. In addition, the simulated TC track largely determined the location of the main rainband. For example, the northern tracks in both the YSU and BouLac schemes and the southern tracks in both the QNSE and MYJ schemes, as mentioned in section 4.1, corresponded to the main rainfall, mostly located more northward in the SR in the former and more southward in the latter. Meanwhile, heavy rainfall occurred in regions with significant low-level convergence, which is known to be an important factor for maintaining precipitation (Zhu et al., 2000). There was a convergent zone between the easterly and northeasterly in the SR in each simulation at 850 hPa (Fig. 11), with the exact location varying with the PBL scheme used. In the SR, the convergences in group 1 were mostly stronger than those in group 2 (figure omitted). Therefore, the difference in the horizontal distribution in precipitation was partly related to the difference in the simulated storm track and partly associated with the low-level convergence distribution.

Figure 10.  Simulated hourly accumulated precipitation (units: mm) from 0600 to 0900 UTC 7 October 2013, using the (a) YSU, (b) MYNN2, (c) BouLac, (d) MYJ, and (e) QNSE PBL schemes. The dashed rectangle (30^°-31^°N, 119.2^°-121^°E) in each panel marks the SR.

Figure 11.  Winds (arrows; units: m s^-1) and divergence field (shaded; units: 10^-5 s^-1) at 850 hPa at hourly intervals for all experiments from 0600 to 0900 UTC 7 October 2013: (a) YSU; (b) MYNN2; (c) BouLac; (d) MYJ; and (e) QNSE. The dashed rectangle (30^°-31^°N, 119.2^°-121^°E) in each panel marks the SR.

Precipitation intensity is closely related to moisture and vertical motion. The hourly vertical profiles of relative humidity averaged in the SR from 0600 to 0900 UTC 7 October 2013, for all simulations, indicated distinct moisture discrepancies in the middle troposphere above 700 hPa (Fig. 12). Group 2 simulated a drier atmosphere than group 1 by as much as 15%. The moist layer with a relative humidity greater than 85% was shallower in group 2 than in group 1 (550 hPa versus 400 hPa at 0900 UTC 7 October 2013). The higher moisture above 700 hPa in group 1 may have been associated with a larger moistening effect of deeper convection due to stronger low-level convergence (Fig. 11), ascending motion (Fig. 13) and vertical mixing (Fig. 14).

Figure 12.  Vertical profiles of the area-averaged relative humidity (units: %) over the SR, marked in Fig. 8, from 0600 to 0900 UTC 7 October 2013, for all simulations.

Figure 13.  Vertical profiles of the area-averaged vertical velocity (units: m s^-1) over the SR region, marked in Fig. 8, from 0600 to 0900 UTC 7 October 2013, for all simulations.

Figure 14.  Vertical-longitude cross sections of the heat diffusion coefficients K_h (shaded; units: m² s^-1) averaged in latitude over the SR region, marked in Fig. 8, at 0700 UTC 7 October 2013 , for all simulations: (a) YSU; (b) MYNN2; (c) BouLac; (d) MYJ; and (e) QNSE.

Moreover, the upward motions were weaker in group 2 than in group 1 (Fig. 13), which was consistent with the stronger convergence in the boundary layer mentioned above. For instance, the maximum upward motion in the YSU scheme (approximately 0.28 m s^-1) was nearly three times that (approximately 0.09 m s^-1) in the MYJ scheme at 0600 UTC 7 October 2013. Again, the vertical motion profiles in group 1, including those in YSU, MYNN2 and BouLac, were quite similar, including the maximum upward motion in the middle troposphere. The major difference among the three schemes lies in the height and strength of ascending motion. The ascending motion was strongest but occurred at the lowest height in the YSU scheme, followed by the MYNN2 scheme. The BouLac scheme produced the weakest ascending motion, except at 0600 UTC 7 October 2013, which was consistent with its weakest simulated rainfall intensity. In addition, from 0700 to 0800 UTC 7 October 2013, the higher altitude of peak vertical velocity for the BouLac scheme may have been related to the deeper neutral and unstable stratification, which allowed the updraft to be accelerated and keep its maximum up to a higher altitude. The true reason may be more complicated and needs to be investigated further.

Additional analyses of the mean surface fluxes over the SR region showed that group 2 produced slightly higher moisture fluxes and less sensible heat fluxes than group 1, resulting in a slightly moister (as indicated in Fig. 12) and colder low-level boundary layer in group 2. However, group 2 led to a significantly weaker averaged vertical eddy diffusivity of heat in the SR than group 1 (see Fig. 14). In the areal mean profile over the SR, the greatest maximum vertical mixing in the boundary layer was produced in YSU, followed by BouLac, MYNN2, QNSE and MYJ, in that order. Therefore, relative to group 2, group 1 featured a more efficient transport of moisture and energy and stronger upward motions in the upper-level PBL and free atmosphere, as indicated in Fig. 13, which were more favorable for the development of a moister atmosphere, convection and precipitation.

The boundary layer structure and rainfall after the landfall of Typhoon Fitow (2013), which caused severe disasters in Zhejiang Province, China, were studied using the Advanced Research version of the WRF model. The study mainly focused on the sensitivity of the simulation to the choice of PBL parameterization scheme in the model. Five PBL schemes (YSU, MYNN2, BouLac, MYJ, and QNSE), separated into group 1 (YSU, MYNN2 and BouLac) with the same surface layer scheme, and group 2 (MYJ and QNSE) with different surface layer schemes, in the WRF model, were examined. The results showed that the simulation with the YSU scheme, called the control run, successfully reproduced the main features of the track, intensity and rainfall of Fitow, with average track and intensity errors of only 19 km and 2.8 m s^-1, respectively. The 6-h accumulated area-averaged rainfall errors were generally less than 10% of the total rainfall, and the spatial correlation coefficients between the simulation and the observation varied from 0.49 to 0.8.

The simulated storm intensity was insensitive to the PBL scheme in the Fitow case, with the mean deviation of the maximum sustained near-surface wind speed among the five simulations less than 3.5 m s^-1. This result is different from those in several previous studies, mainly due to the quick weakening and short time integration for the landfalling TC in this study. However, the simulated storm track after landfall was sensitive to the PBL scheme, with a deviation as large as 115 km at one point. This deviation was found to be related to the difference in the simulated steering flow among the PBL schemes.

The simulated boundary layer structure was very sensitive to the choice of PBL and surface layer scheme. The difference in the average maximum tangential (radial) wind among all simulations was as large as 11% (33%) of that in the simulation with the YSU scheme. Group 1 simulated stronger vertical mixing than group 2. This result explains the simulated weaker vertical shears of tangential and radial winds as well as the deeper inflow layer and the well-mixed boundary layer in group 1.

The simulated rainfall distribution associated with Typhoon Fitow (2013) varied greatly among different PBL schemes, especially after landfall. The maximum hourly rainfall was as much as three times the minimum in the SR between group 1 and group 2. The differences were related partly to the difference in the simulated storm track and partly to the difference in the simulated low-level convergence and vertical motion among different PBL schemes. Group 1 produced stronger vertical mixing, a deeper moist layer and stronger ascending motion than group 2, which were beneficial for stronger convection and larger rainfall development in the SR.

An interesting finding is that the differences in both boundary layer structure and rainfall among simulations with different surface layer schemes were considerably larger than those among simulations with the same surface layer scheme. This finding suggests that a more accurate representation of surface layer fluxes is key to the realistic simulation and prediction of the boundary layer structure and rainfall associated with landfalling TCs.

Although our results show better performance in the simulation of Typhoon Fitow (2013) with the YSU scheme, this result may not be broadly conclusive because we conducted only one case simulation. In addition, because of the nonlinearity of physical processes and their feedbacks, the sensitivity could vary with different combinations of other model physics parameterizations. Nevertheless, our results strongly suggest that caution needs to be taken when one chooses a surface layer and PBL scheme for the simulation and prediction of landfalling TCs. Realistically simulating/predicting heavy rainfall associated with landfalling TCs remains a great challenge. More studies should be devoted to the improved understanding and prediction of severe weather associated with landfalling TCs. In particular, it is important in future studies to further identify the strengths and weaknesses of the available physical parameterization schemes and to develop new parameterization schemes suitable for skillful TC simulation/prediction.