A Modified Double-Moment Bulk Microphysics Scheme Geared toward the East Asian Monsoon Region

Cloud processes and their associated feedbacks play significant roles in weather and climate models; thus, cloud microphysical parameterization development has been a major research field for the last several decades. The Intergovernmental Panel on Climate Change (IPCC) confirmed that cloud-related processes remain the largest contribution to the overall uncertainty in the present numerical models (IPCC, 2021). One of the main reasons for the uncertainty is that numerous microphysical and dynamical processes controlling the life cycle and radiative properties of clouds are not correctly represented in the numerical models.

Combining the effects of the Tibetan plateau, East Asia Monsoon, North Pacific Ocean, and intensive human activities, the physical characteristics of cloud and precipitation over the East Asian monsoon region (EAMR) are not typically shared by other monsoon systems (e.g., Wang et al., 2004; Xu, 2012; Yin et al., 2013a; Li et al., 2016; Chen et al., 2017), which has resulted in troublesome issues in terms of weather prediction and climate projection in this region (e.g., Jiang et al., 2015; Park et al., 2016; Kim et al., 2018; Yao et al., 2018). So far, more than 30 cloud microphysics schemes with varying degrees of sophistication, from the simple warm scheme to the complicated bin microphysical scheme, have been developed. Because of limited observations on cloud microphysics over the EAMR, most existing microphysics schemes were developed based on observations outside East Asia. Thus, the existing microphysics schemes did not consider cloud and precipitation processes unique to the EAMR. For example, the high aerosol concentration over the EAMR, which may play significant roles in clouds and precipitation (e.g., Qian et al., 2009; Fan et al., 2012; Guo et al., 2017; Lin et al., 2018), has been not considered in the schemes. One of the possible effective ways to improve the numerical prediction of clouds and precipitation over the EAMR is to optimize the microphysical schemes based on the observations from this region.

With the intent of developing a microphysics scheme geared toward the unique cloud and precipitation physics over EAMR, we have analyzed cloud and precipitation microphysical properties in this region over the last decade using multi-source observational datasets, including long-term in-situ measurement and remote sensing datasets (Yin et al., 2011, 2012, 2013a, b; Yin et al., 2014a, b; Wang et al., 2015; Yin et al., 2015; Wang et al., 2019). Based on the observations and theoretical research (e.g., Xu and Wang, 1985; Xu and Duan, 1999; Xu and Yin, 2017; Yin et al., 2022), an improved Milbrandt-Yau (named IMY) scheme has been developed based on the standard Milbrandt-Yau (MY hereafter) scheme in the Weather Research and Forecasting (WRF) model version 3.5.

The remainder of the paper is structured as follows. General properties of the MY scheme are presented in section 2. Section 3 focuses on the detailed improvements of cloud microphysical processes. A description of two heavy rainfall cases and model configuration is presented in section 4. Section 5 focuses on the simulated results. Conclusions and a summary are given in section 6.

To the best of our knowledge, the MY scheme (Milbrandt and Yau, 2005a, b) is a relatively complete scheme with various cloud microphysical processes, and it is easy to add or modify a cloud-precipitation microphysical process. Consequently, we adopt the MY scheme as a foundation and develop a series of modifications. Interested readers are requested to refer to the responding references for a more detailed description of each process (Ferrier, 1994; Milbrandt and Yau, 2005a, b). Overall, two liquids (cloud droplet and raindrop) and four ice hydrometeors (cloud ice, snow, graupel, and hail) are included in the MY scheme. According to the definition, liquid water particles with a diameter (D) larger than 200 μm are classified as raindrops, while those having a D less than 200 μm are defined as cloud droplets. Cloud ice is a pristine ice crystal with a maximum dimension length of less than 300 μm, and snow contains large ice crystals with a length greater than 300 μm. Graupel (hail) is a moderate (high) density solid particle. For simplicity, graupel (hail) is usually limited to D or its maximum dimension length less (greater) than 5 mm. The densities of cloud and rainwater (${\rho _{\text{w}}} $), ice (${\rho _{\text{i}}} $), snow (${\rho _{\text{s}}} $), graupel (${\rho _{\text{g}}} $), and hail (${\rho _{\text{h}}} $) are 1000 kg m^–3, 500 kg m^–3, 100 kg m^–3, 400 kg m^–3, and 900 kg m^–3, respectively.

It is assumed that the mass (m) of a particle is related to its D by

where x refers to cloud water, rain, cloud ice, snow, graupel, or hail. Detailed information of constant parameters c_x and d_x are listed in Table 1. It should be noted that ice and snow particles are not always spherical owing to their complicated structure. Therefore, the observation-based fitted parameters in Field et al. (2005) and Brandes et al. (2007) are employed for cloud ice and snow, respectively, in the IMY scheme. Notably, options that assume that ice and snow are spherical are retained.

All hydrometeor categories have nonzero fall speeds except for cloud droplets. The relationship that represents the terminal fall velocities for each category is given by

Here, ${\rho _0} $ is the reference air density (usually at 1000 hPa level), $\rho $ is the air density, and the air density correction factor is given by ${\left( {{\raise0.7ex\hbox{${{\rho _0}}$} \mathord{\left/ {\vphantom {{{\rho _0}} \rho }}\right.} \lower0.7ex\hbox{$\rho $}}} \right)^{0.5}} $. The values of the coefficients a_x, b_x, and f_x for each hydrometeor species are listed in Table 2. It is worth noting that the coefficients significantly influence cloud microphysical processes and surface rainfall intensity. In light of this, different coefficients are introduced as options in the IMY scheme, which offers users the flexibility to choose an option expediently. According to our previous case study (Zhang et al., 2016), the parameters from Ferrier (1994) are suitable for cloud and precipitation microphysics over north China, which is set as the default option.

The gamma size distribution is utilized to represent particle size distribution for each hydrometeor category in the MY scheme, which is expressed by

Here, N_x (D_x) is the concentration per unit volume of particles of diameter D_x and D_x + δD for category x, ${\alpha _x} $ is shape parameter, and ${\lambda _x} $ is interpreted as slope. N_0x is often referred to as the intercept and shows a wide range of variation with dependence to ${\lambda _x} $ and ${\alpha _x} $, and it is given by

Here, N_Tx is the total number concentration. Combining particle size distribution [Eq. (3)] and the relationship between mass and diameter [Eq. (1)], the slope of a given hydrometeor category is obtained by

Here, q_x denotes the mixing ratio. Incidentally, the exponential size distribution, widely used to represent hydrometeor particle size distribution in the early days, can be obtained as a particular case of the gamma size distribution with ${\alpha _x} $ = 0. The following section highlights the shape parameters for cloud droplet, raindrop, and snow size distributions.

The shape parameter plays an important role in the evolution of particle size distribution and cloud microphysical processes (Zhang et al., 2003; Milbrandt and Yau, 2005b; Seifert, 2005; Tao and Hong, 2009; Milbrandt and McTaggart-Cowan, 2010; Luo et al., 2021; Milbrandt et al., 2021). In the MY scheme, the shape parameters for cloud droplet, raindrop, and snow particle size distributions are fixed values of 3, 0, and 0, respectively. In reality, shape parameters vary across a wide range in nature, and thus fixed values would cause large limitations to particle size distributions. In the IMY scheme, however, those values are permitted to vary based on the previous observations rather than existing as fixed constants.

According to long-term observations (Yin et al., 2011), the shape parameter (${\alpha _{\text{c}}} $) within a cloud droplet size distribution spans a wide range of 0–12, and thus a fixed value for ${\alpha _{\text{c}}} $ would cause large errors within the cloud droplet size distribution. In the IMY scheme ${\alpha _{\text{c}}} $ is diagnosed from the cloud number concentration (N_c), following Thompson (2008), with a modification that is based on the long-term observations over EAMS, which is given by

Zhang et al. (2003) proposed that the ${\alpha _{\text{r}}} $–${\lambda _{\text{r}}} $ relation can capture the mean physical characteristics of raindrop spectra. Numerous studies (e.g., Brandes et al., 2003; Yang et al., 2010; Chen et al., 2013) provided a relationship between ${\alpha _{\text{r}}} $ and ${\lambda _{\text{r}}} $ for a raindrop size distribution, which is a useful method to retrieve ${\alpha _{\text{r}}} $ according to the prognostic value of ${\lambda _{\text{r}}} $. In the IMY scheme, we retrieve the shape parameter by employing the empirical relation between ${\alpha _{\text{r}}} $ and ${\lambda _{\text{r}}} $ based on the statistical results from numerous measurements in the past several decades over the EAMR (Wang et al., 2015), with the expression being given by

Similar to the shape parameter of raindrop size distribution, potential relationships between slope and shape for the snow size distribution have been proposed in the last decade (e.g., Zhao et al., 2010). In the IMY scheme, the shape parameter for snow size distribution is derived from the empirical relation between ${\alpha _{\text{s}}} $ and ${\lambda _{\text{s}}} $ given by Zhao et al. (2010), which is expressed as

It should be noted, like any constant parameter in the microphysics schemes, the shape-slope relationship obtained over the Qilian Mountains may not well represent all the conditions over the East Asian monsoon region. Users may adjust the shape-slope relationship according to upcoming observations if this approach works.

Great care should be taken to apply a choice of terminal fall velocities because it directly affects the vertical flux of hydrometeors and even local air motions. There are several ways to represent the terminal fall velocities used in the collection equations and sedimentation, including the mass-weighted mean, number-weighted mean, diameter-weighted mean, and mean velocity. Usually, the mass-weighted mean fall velocity (V_q) is widely utilized for the hydrometeor mixing ratio in a single moment microphysics scheme. It is worth noting that Xu and Duan (1999) proposed that the mass-weighted mean terminal velocity is much larger compared to the mean terminal velocity.

A second predictive equation for number concentration (N_x) is included in a double-moment scheme, and number concentration-weighted fall velocity (V_N) is also introduced. Milbrandt and Yau (2005a) provided an in-depth discussion on the advantages of the double-moment over the single-moment approach. They proposed that size sorting can be resolved by the double-moment schemes that include a second predictive equation for the number-weighted mean terminal velocity. Size sorting, which has strong effects on cloud microphysical processes, is one of the features in a double-moment scheme with different sediment fall velocities of V_q and V_N. One serious problem is the uncertainty of N_x, which varies over a large range for a given q_x. In addition, Milbrandt and Yau (2005a) state that given enough time, size sorting can eventually cause unrealistically large mean sizes because the fall speed ratio always exceeds one. Consequently, special numerical artifacts were introduced to limit the problem (e.g., Cohard and Pinty, 2000; Milbrandt and Yau, 2005a; Milbrandt and McTaggart-Cowan, 2010). To the best of our knowledge, the predicted number concentration has a significantly larger uncertainty than the mixing ratio because droplet number concentration is controlled by complicated factors (Van Weverberg et al., 2012). Therefore, we use the same terminal fall velocity V_q for both q_x and N_x sedimentation in the IMY scheme at present, and V_q is given by

It should be noted that size sorting is misapplied to both the mixing ratio and number concentration of hydrometeors. Keeping this shortcoming in mind, we have launched future work to treat the size sorting in the IMY scheme. To avoid the uncertainty in V_N, we still use V_q but with a turning factor whose value is less than 1.0, as V_N is always less than V_q. A series of tests were launched, and the results were compared to polarimetric radar observations. This approach will be introduced in the IMY scheme as an option in future applications.

Based on field experiments measurements, we investigated ice nuclei properties and compared them with other regions (Yin et al., 2012; 2013a). In the IMY scheme, high ice nuclei concentrations over East Asia are considered, and the expression of ice nucleation is written as

Similar to Reisner et al. (1998), we have limited the value of T in the above equations for temperatures cold than 243 K because the Fletcher-type curve produces too many ice particles at temperatures below 243 K. Cloud ice mixing ratio tendency is obtained with an assumption of the original nucleated ice particle mass of ${m_{{\text{i}}0}} $ = 10⁻¹² kg, giving

In the IMY scheme, the autoconversion process of cloud water to form rainwater has been replaced with an improved Kessler-type parameterization given by Yin et al. (2015). One of the novel features of their parameterization is that the threshold is diagnosed as a function of cloud altitude. The autoconversion from cloud water to rainwater is given as

Here, $\alpha $ = 0.001 s^–1 is a time constant, H represents the Heaviside function, and the threshold ${q_0} $ varies with the altitude of a cloud. The corresponding equation is shown as

Here, h is the altitude of a cloud above sea level, and ln represents the natural logarithmic function. Note that the unit of ${q_{\text{0}}} $ is g m^–3.

In the MY scheme, the Hallet-Mossop processes (Hallet and Mossop, 1974) were utilized to represent ice multiplication, which only occurs under the limited temperature range (265–270 K) and in the presence of supercooled water droplets larger than 24 μm. To expand the Hallet-Mossop processes in the IMY scheme, we introduce ice multiplication due to ice–ice collision processes, according to Yano and Phillips (2010). The collisions among hail, graupel, and snow were considered, and the corresponding equations can be expressed as follows:

Here, ${\alpha _0} = 1.2 \times {10^{ - 3}}\;{{\text{m}}^3}\;{{\text{s}}^{ - 1}} $ is a rate coefficient. The increasing mass of cloud ice associated with this process is given as

It is necessary to retain a saturation adjustment technique in non-hydrostatic models, which calculates the amount of condensation and/or deposition to remove any supersaturated moisture, or the amount of evaporation and/or sublimation to remove any sub-saturation in the presence of cloud droplets and cloud ice (Tao et al., 1989). In the MY scheme, a simple saturation adjustment proposed by Kong and Yau (1997) was used, which only took liquid water into account. In the IMY scheme, an ice-water saturation adjustment scheme given by Tao et al. (1989) is introduced; in this scheme, there are two major assumptions, (1) the saturation vapor mixing ratio is defined as a mass-weighted combination of the saturation values over liquid water and ice within the temperature range of the triple point of water (T₀) and the temperature of homogeneous freezing (T₀₀); and (2) condensation and deposition occur in proportions that depend linearly on the temperature in the range T₀ to T₀₀. Water vapor would be assumed to condense/deposit into cloud water/ice under a supersaturated condition. On the contrary, cloud water/ice would be assumed to evaporate/sublimate when subsaturation occurs. Note that large drops (raindrop, snow, and graupel) evaporate and/or sublimate only in cases of subsaturation after all cloud water and cloud ice are exhausted.

In a double-moment scheme, q_x and N_x are predicted independently. As a result, N_x can vary over a considerable range (up to several orders of magnitude) even with the same q_x. Such an approach may lead to a negative or very large value of N_x for a given q_x, and thus would cause an imbalance between mass and number and even erroneous results (Yin et al., 2014b). To the best of our knowledge, most of the existing double-moment microphysics schemes do not include any constraint condition to limit N_x within tolerable bounds. Sometimes, microphysical variables may attain negative values, and when such cases did occur, the negative values were always simply set to zero (e.g., Reisner et al., 1998). Additionally, an upper bound value is given for convenience in the microphysical parameterizations since the independent change of q_x on N_x is implied (Morrison and Pinto, 2005; Gao et al., 2011).

We have noticed that N_c was shown to attain a negative value in several real case simulations (e.g., Wang et al., 2014). Consequently, a relationship between q_c and N_c has been proposed within a context adhering to the necessity of the cloud microphysical parameterization (Yin et al., 2014b). The relation between q_c and N_c is given by

Another issue centers around the assumption that all cloud ice particles immediately melt at the level where the temperature is greater than 0°C. Thus the latent heating from cloud ice particle melting sometimes leads to obvious oscillations around the melting point and associated temperature levels. Given this shortcoming, numerical artifacts were introduced with the latent heating redistributed in the vertical levels. In this way, 60% of the total latent heating is released at the melting level, 30% at the first lower level, and the remainder at the third level (downward). If the melting level is the second model level, 70% of the total latent heating is released at the melting level, and 30% is released at the lower one. It should be noted that the numerical artifact approach was only performed for the latent heating from cloud ice melting. Large ice particles (i.e., snow, graupel, and hail) melt gradually in thermodynamic equilibrium with the surrounding air.

A simulation based on real data is a robust approach to understand the performance of a cloud microphysics scheme by comparing modeled products against observations, such as surface rain gauge, radar, and satellite observations. Two heavy rainfall cases were selected to validate the IMY scheme in the study. One was a mid-latitude extreme rainfall event in northern China, and the other was a local extreme rainfall event in the warm sector over southern China.

The heaviest rainfall in Beijing since 1951 occurred on 21 July 2012, with a maximum rainfall of 460 mm recorded at the Fangshan gauge station, Beijing, China. Several rain gauge stations recorded hourly rainfall exceeding 50 mm, and the average rainfall throughout Beijing City was 170 mm (Zhao et al., 2013). Although most of the current operational models were able to reasonably predict the extreme rainfall event one day prior, there were notable timing and location errors, and the amount of the rainfall was significantly underestimated by the operational predictions, especially the rainfall in the warm sector ahead of the cold front (Zhang et al., 2013).

Figure 1 displays a weather chart at 0000 UTC on 21 July 2012, 2–3 h before the rainfall event. The meridional circulation at the 500 hPa level was characterized by two blocking highs over the Ural Mountains and the Okhotsk Sea. Previous studies (e.g., Zhong et al., 2015; Mao et al., 2018) have confirmed that the heavy rainfall event was closely related to the low-pressure center situated over the Baikal Lake. From the 850 hPa level, one can see that a branch of water vapor was transported to northern China by the southwest and southeast airflows, which was closely related to the typhoon Vicente over the South China Sea.

Figure 1.  Composite weather charts at 0000 UTC 21 July 2012: geopotential height (black-contoured at 20 m intervals) and temperature (red contoured at 4°C intervals) at 500 hPa, and wind barbs (a full barb is 4 m s^–1) at 850 hPa, and precipitable water (kg m^–2, shading). Nested model domains used for the WRF simulations are marked by black rectangles with grid resolutions of 36 km (D01), 12 km (D02), and 4 km (D03), respectively. The circle-cross marker denotes the location of the radiosonde site in Beijing.

The WRF model was configured with three nested domains of 36 km, 12 km, and 4 km horizontal resolution. The ERA-Interim reanalysis datasets (Dee et al., 2011) with 0.75-degree spatial resolution at temporal intervals of six hours were adopted for initialization and boundary conditions. To improve the initialization conditions, conventional observation data, including sounding and surface observations, were assimilated by the Gridpoint Statistical Interpolation (GSI) data assimilation system (Kleist et al., 2009). The geographical coverage of each domain is shown in Fig. 2, and a total of 57 sigma levels were assigned vertically with the model top up to 20 hPa. There are 15 levels below 850 hPa and 9 levels above 200 hPa. The model was initiated at 0000 UTC on 21 July 2012 and ran for 24 hours, and two-way feedback is applied in configuration. Observational rainfall was taken from the hourly automatic weather station observation, which contains more than 50 000 stations over China. Ground-based radar observations around Beijing were utilized to generate composite radar reflectivity. For the cloud coverage, the temperature of blackbody brightness (TBB) products from the Chinese geostationary meteorological satellite Fengyun-2E (FY-2E) satellite were used to denote the spatial distribution and intensity of a convective system.

Figure 2.  Same as Fig. 1 but for the extreme local rainfall case in the warm sector over southern China at 0000 UTC 8 May 2014. The circle-cross marker indicates the location of the radar station in Guangzhou, China.

General information of the model configuration is summarized as follows. The Kain-Fritsch (new Eta) cumulus parameterization scheme (Kain, 2004) has been utilized for the outer two domains, while no cumulus scheme is used for the finest domain. The rapid radiative transfer model (rrtm) scheme (Mlawer et al., 1997), the Mellor–Yamada–Janjic Scheme (MYJ) scheme (Janjić, 1994), the revised MM5 Monin-Obukhov scheme (Jiménez et al., 2011), and the unified Noah land surface scheme (Tewari et al., 2004) have been applied for all domains.

An extreme local rainfall event occurred in the warm sector over southern China on 8 May 2014, which caused casualties and property losses. More than ten rain gauge stations recorded 24-h accumulated rainfall exceeding 100 mm, and a maximum value of 180 mm occurred at Shanwei station, Guangdong Province, China. The heavy rainfall belt was mostly along with the coastline of southern China. The heavy rainfall occurred in the warm sector, far from the cold frontal zone. Operational forecasts almost missed the heavy rainfall at that time. Figure 2 displays a weather chart at 0000 UTC 8 May 2014. Southern China was situated on the western edge of the subtropical high marked by a 5880 gpm isohypse at 500 hPa. Meanwhile, two low-height centers were obvious over northern China. However, the low centers were far away from southern China, and thus they had no direct effect on the rainfall event due to the long distance. It is noteworthy that southern China was dominated by strong southwesterly or southeasterly flow, which was favorable for the heavy rainfall because of the large expanse of warm moisture-laden air transported into southern China. Nevertheless, large amounts of precipitable water (>60 kg m^–2) were present over southern China and the South China Sea. Owing to the large precipitable water core near Hainan Island, China, warm moist air was advected into Guangdong by the strong southwesterly flow.

The model configuration is similar to the mid-latitude extreme rainfall case with slight modifications. The geographical coverage of each domain (Fig. 2) is as follows. The 4 km domain (D03) covers southern China with grid points of 505 × 274, surrounded by a 12 km domain (D02) with 373 × 232 grids. The 4 km and 12 km domains are nested within a 36 km domain (D01) with grid points of 270 × 192. It should be noted that the ERA-interim reanalysis data have a large bias of water vapor over southern China, which had a significant effect on the simulations. In light of this shortcoming, the NCEP final analysis (NCEP-FNL), consisting of one-degree datasets at intervals of six hours, was utilized as the model initialization and boundary conditions. Owing to the differences in surface properties between northern China and southern China, the Yonsei University (YSU) scheme (Hong et al., 2006) was applied for the boundary layer and the Noah-MP scheme (Niu et al., 2011) was introduced for land surface processes. The WRF model was integrated for 24 h, starting at 0000 UTC 8 May 2014, with outputs at 10-min intervals.

Observed and simulated 24-h accumulated rainfalls are shown in Fig. 3. Generally, the event was well-simulated using the WRF model with the IMY scheme. More specifically, the model successfully replicated the entire rainfall band distribution, and the simulated rainfall bears a substantial resemblance to observation. Specifically, the model produced heavy rainfall of over 200 mm in southwestern Beijing. Notably, a local area of strong precipitation located northeast of Beijing was well replicated. Although the main properties of the rainfall were well reproduced, it appears that the MY scheme significantly underestimated the heavy rainfall core.

Figure 3.  Observed (black contours) and simulated (shadings) 24-h accumulated rainfall (mm) during the period from 0000 UTC 21 to 0000 UTC 22 July 2012, in which (a) and (b) denote the simulations with the IMY and MY schemes, respectively. The rectangles (39.5°–40.3°N, 115.3°–116.2°E) marked with white lines are used for domain averaging in Fig. 5.

Figure 4 illustrates the root mean square error (RMSE) of the 24-h accumulated rainfall by using the surface rain gauge observations with 3141 stations in total located in Fig. 3. Comparatively speaking, the IMY scheme has less RMSE than the MY scheme. It is worth noting that the RMSE of the rainfall over 250 mm is significantly reduced. Figure 3 shows that the MY scheme missed most rainfall over 250 mm.

Figure 4.  The RMSE of 24-h accumulated rainfall (mm) from the simulations with the IMY and MY schemes from 0000 UTC 21 to 0000 UTC 22 July 2012 (See text for details).

To reveal features of the timing variability of the simulation rainfall for this heavy rainfall event, both the observed and simulated area-averaged hourly rainfall rates over the mainly heavy rainfall region are shown in Fig. 5. From the observations, one can see that the observed hourly rainfall rates gradually increased from 0100 UTC, reaching their first peak at 0600 UTC with a value of 20.8 mm. After a slight decrease, the hourly rainfall rate reached its second peak of 21.4 mm at 1100 UTC. Then, the hourly rainfall rates gradually decreased and ceased at 2123 UTC. As for the simulations with the IMY scheme, the temporal evolutions of the simulated hourly rainfall rates from the nested domains are close to the observations, appearing as a double peak mode. The maximum simulated peak hourly rainfall rate of 23.1 mm h^–1 from the former is slightly larger than the latter of 20.4 mm h^–1, and both are close to the observed value of 21.4 mm h^–1. However, the MY scheme fails to replicate the observed double peak mode and underestimates the hourly rainfall rates.

Figure 5.  Time series of domain-averaged hourly rainfall rates (mm) from the different grid spacing domains and rain gauge observations during the period from 0000 UTC 21 July 2012 to 0000 UTC 22 July 2012, in which (a) and (b) denote the simulations with IMY and MY schemes, respectively. The locations of the averaged areas are enclosed by the white lines in Fig. 3.

Observed and simulated composite radar reflectivity (CR) at 0530, 1100, and 1830 UTC on 21 July are presented in Fig. 6. From the observations, one can see that strong radar reflectivity first appeared at 0530 UTC in the southern area of Beijing, which corresponded to the first strong rainfall stage. At 1100 UTC, a strong radar reflectivity band (>40 dBZ) crossed Beijing from the northeast to the southwest, and the main convective band moved toward the east and out of Beijing after 1830 UTC. As for the simulated CR, the model reproduced acceptable radar reflectivity structures and their evolutions. More specifically, the simulations successfully captured the features related to the rainfall event's formation, development, and disappearance. However, the simulated radar reflectivity over 40 dBZ was slightly more extensive than what was observed. One of the potential sources for the difference might be caused by radar observation data processes. Some information is over-treated when radar data from several stations are combined with strict quality control (Wang et al., 2009). Although there are some differences, both the IMY and MY schemes can capture the main features of the convective system. For example, the former strong radar reflectivity belt at 1100 UTC is closer to the observed, compared to that of the latter (Figs. 6b, e, h). In addition, obvious false radar reflectivity over the southwestern Beijing was generated by the MY scheme (Fig. 6i).

Figure 6.  Horizontal maps of composite radar reflectivity from (a–c) observations, and (d–f) and (g–i) simulations with IMY and MY schemes, respectively, at 0530, 1100, and 1830 UTC 21 July 2012.

Now we move to a zoomed-out view of the whole rain belt by comparing the spatial distribution of simulated TBB to the observed values at 0500, 1100, and 1800 UTC 21 July 2012 (Fig. 7). Observations revealed a convective band formed at approximately 0500 UTC extending from the northeast to southwest China. As the convective band moved eastward, several organized convective cells were generated, characterized by TBBs colder than –70°C, noting that the Unified Post Processor (UPP) from the National Centers for Environmental Prediction (NCEP) was applied to retrieve the simulated TBB. For example, four organized convective systems were obvious at 1800 UTC. The convective band and its evolution were successfully reproduced in both simulations. One can see that the IMY scheme successfully reproduced several embedded convective cells. In contrast, the organized convective cells were indiscernible in the simulation with the MY scheme (Figs. 7c, f, i). Although the convective band was duplicated successfully, the simulated coverage is somewhat larger than the observation. Similarly, the retrieval processes and model horizontal resolution might have led to some differences between the simulation and the observation.

Figure 7.  Horizontal maps of the temperature of blackbody brightness (TBB) from (a–c) observations and (d–f) and (g–i) simulations with IMY and MY schemes, respectively, at 0500, 1100, and 1800 UTC 21 July 2012.

Upon comparing simulated radar reflectivity and TBB to those observed, the results show that the convective band and embedded cells and their evolutions were reproduced successfully. The simulations with the IMY scheme are comparable to the observations, while the MY scheme seems to underestimate the convective band and does not resolve the convective cells within the convective band.

Time-height cross-sections of area-averaged hydrometeors (of cloud water, rain, cloud ice, snow, graupel, and hail) mixing ratios are shown in Fig. 8. For the simulation with the IMY scheme, it is apparent that there was large cloud water content during the heavy rainfall process, with a peak value of 0.45 g kg^–1 during the period of 0600 UTC to 0800 UTC. Cloud water concentrated between the 500 hPa and 300 hPa levels. Rainwater was mainly located below 600 hPa with the maximum rain mixing ratio of 0.7 g kg^–1, and the rainwater reached the ground. Cloud ice mainly appeared above the height of 300 hPa, with a maximum value greater than 0.3 g kg^–1. The spatial distribution of snow was close to cloud ice, while the snow content was far less than that of cloud ice. Graupel mainly occurred above the height of 600 hPa, and most of it was concentrated between 500 and 300 hPa levels, with the maximum value of graupel mixing ratio greater than 0.45 g kg^–1. There was a very low hail mixing ratio, similar to snow, and the hail mixing ratios were less than 2.1 × 10^–3 g kg^–1 and were mainly located below 500 hPa. In summary, there were large mixing ratios of cloud water, rain, cloud ice, and graupel, while the mixing ratios of snow and hail were small. Cloud ice, snow, graupel occurred at the upper levels, and rain and hail were mainly concentrated at lower levels. Cloud water had a wide spatial distribution. These results helped outline three-dimensional cloud structures of the heavy rainstorms over northern China.

Figure 8.  Time-height cross-sections of the area-averaged hydrometeor mixing ratios from simulations with (a–g) the IMY scheme and (h–n) the MY scheme. The averaging area (39.0°–41.0°N, 114.0°–117.0°E) nearly covers the entire area with rainfall greater than 100 mm.

Hydrometeor cross-section patterns produced by the MY scheme are similar to those given by the IMY scheme. However, differences are visible between them except for water vapor. The most obvious differences are the cloud ice and snow, followed by graupel and rainwater, with a little difference regarding cloud water (Figs. 8b–h). The increased cloud ice may correspond to the high ice nuclei concentration over northern China. Under a rich water vapor supply condition, a higher ice nuclei concentration would lead to higher numbers of cloud ice particles. Sequentially, the higher cloud ice is favorable for generating large ice particles, including snow, graupel, and hail, by cold cloud microphysical process. Consequently, the enhanced cold cloud processes contribute partially to the rainwater (Fig. 8c) and the surface rainfall (Fig. 3).

Modifications in the cloud microphysical scheme also change the dynamical fields because cloud microphysical and dynamical processes are entangled in the clouds. The changes to the updraft can be seen from a cumulative contoured frequency by altitude diagram (CCFAD) given in Fig. 9. The CCFAD presents the percentage of horizontal grid points with vertical motion weaker than the abscissa-scaled value for a given height (Yuter and Houze, 1995). Here, vertical speeds are binned with intervals of 1 m s⁻¹ based on 13 model outputs with half-hour intervals during the severe rainfall episode from 0600 UTC to 1200 UTC 21 July 2012. Overall, the IMY scheme provides similar CCFAD patterns of updraft motion below 7 km compared to those given by the MY scheme, while differences can be seen in the upper levels. One can see that the IMY scheme generates stronger updraft motions above 7 km than those of the MY scheme. This is why the IMY scheme shows stronger ice cloud microphysical processes in the upper levels than those of the MY scheme (Fig. 8).

Figure 9.  The cumulative contoured frequency by altitude diagram (CCFAD) of the simulated updraft motion for the IMY scheme (shadings) and the MY scheme (contours) within the respective boxes marked white lines in Fig. 3. The CCFAD is calculated from thirteen model outputs with half-hour intervals during the severe rainfall episode from 0600 UTC to 1200 UTC 21 July 2012.

Comparisons between observed and simulated 24-h accumulated rainfall are illustrated in Fig. 10. From the observations, one can see that the heavy rainfall belt with values over 100 mm was distributed along the coastline of Guangdong Province. The rainfall amounts decreased away from the coastline toward the northwest, while several rainfall centers were located near 26°N with values over 25 mm. Overall, the spatial distribution of heavy precipitation was well reproduced. The simulated heavy rainfall belt along the coastline shows good agreement with observations. The local rainfall near 26°N was also reproduced, although the rainfall was slightly overestimated. The model simulated a peak rainfall of 204 mm, close to the maximum observational value of 180 mm at the Shanwei station. Specifically, the MY scheme overestimated the rainfall, although the main rainfall belt along the coastal line was well replicated.

Figure 10.  Observed (black contours) and simulated (shadings) 24-h accumulated rainfall (mm) during the period from 0000 UTC 8 to 0000 UTC 9 May 2014, in which (a) and (b) denote the simulations with the IMY and MY schemes, respectively. The rectangles marked with red lines denote averaging areas in Figs. 14 and 15.

Figure 11 shows the root mean square error (RMSE) of the 24-h accumulated rainfall by using the surface rain gauge observations (2593 stations in total) located within the domain of Fig. 10. One can see that the IMY scheme has less RMSE than the MY scheme, except for rainfalls between 10 mm and 25 mm. Compared to the mid-latitude extreme case, smaller improvements are obtained in the local extreme rainfall event in the warm sector over southern China.

Figure 11.  Same as Fig. 4 but for the local extreme rainfall event in the warm sector over southern China.

Observed and simulated composite reflectivities (CRs) are illustrated in Fig. 12. The observations show three obvious CR belts at 0600 UTC from the South China Sea inland to Guangdong Province, China. More specifically, one belt was situated over southern China along the coastal line from the southwest toward to northeast. The other two belts were located in a parallel fashion, from west to east, over the South China Sea. Almost four hours later, at 0948 UTC, a convective band formed and moved slowly toward the southeast. The slowly-moving convective band lasted several hours and produced heavy rainfall along the coastline. The spatial patterns of the simulated CR agree well with the observations. However, the simulated convective band in the middle and later periods (shortly after 0600 UTC) developed slightly faster than observation. As a result, there were some differences in the simulated location and timing. To sum up, the WRF model with the IMY scheme can simulate the development of a major convective system, evidenced by a comparison of the model output with observations from ground-based rain gauges and weather radars.

Figure 12.  Horizontal maps of composite radar reflectivity from (a–c) observations and (d–f) and (g–i) simulations with IMY and MY schemes, respectively. Black lines in panels (b, e, h) denote the locations of vertical cross-sections used in Fig. 13, and the star indicates the location of the Heshan radar station, Guangdong Province, China. Note that the simulated convective band in the middle and later periods (nearly after 0600 UTC) developed slightly faster than the observed. Thus there were some timing differences between the observations and simulations.

Next, we focus on the unique vertical structures through a convective cell by comparing simulated radar reflectivity vertical profiles to C-Pol radar RHI scans observed at a 236° azimuth from the Heshan station on 0952 UTC 8 May 2014 (refer to the radar station and the locations of the vertical cross-section in Fig. 12b). It can be seen that the C-Pol radar RHI scans nearly crossed a convective cell, and the fine internal structures were captured. Overall, the convective cloud (strong radar reflectivity) was mainly located below 5 km, indicating that warm rain processes were dominant. Besides, the convective cell was surrounded by stratiform clouds. Given the vertical profile (Fig. 13), warm rain processes contributed the most to the precipitation. Generally speaking, these special vertical structures were well reproduced by the IMY scheme, although there was a spatial shift. However, the MY scheme overestimated the ice cloud processes. It should be noted that there was a spatial shift between the simulation and observation due to the simulated position errors of the convective cells. As a result, the particular spatial distribution of the MY simulated convective core is not completely consistent with observation.

Figure 13.  Radar reflectivity vertical profiles derived from (a) the C-Pol radar RHI observation at 0952 UTC 8 May 2014, and (b) and (c) simulations with the IMY and MY schemes at 0950 UTC 8 May 2014, respectively. Blue lines in panels (b) and (c) indicate the 0°C isotherm.

Vertical profiles of domain-averaged mixing ratios of water vapor (q_v), liquid water content (q_c plus q_r), and solid hydrometeors (the sum of q_i, q_s, q_g, and q_h) are illustrated in Fig. 14. One can see that there was plenty of water vapor during the whole rainfall process with only very slight variations (Fig. 14a). The amount of water vapor decreased from the surface to 8 km, with the maximum value of water vapor greater than 16 g kg^–1 (Fig. 14d). There were two cores with large liquid water content. The lower core was mainly below 600 hPa (nearly 4 km) with the largest value of 0.72 g kg^–1. The melting of ice particles made a considerable contribution to the total liquid water. Furthermore, rainwater decreased from 700 hPa downward to the ground because of the evaporation process (Fig. 14e). Meanwhile, the upper core was situated between 500 hPa and 300 hPa with a peak hydrometeor mixing ratio of over 0.45 g kg^–1, mostly composed of cloud water. The small liquid drops were transported to upper levels owning to the strong updraft, and thus they could not completely freeze in a short time (Xu et al., 2011). The MY scheme produced similar vertical patterns for water vapor (Fig. 15a); however, it did not generate a liquid water core in the upper levels (Fig. 15b). The IMY scheme permitted considerable liquid water (cloud water) content in levels above 500 hPa (Fig. 14e). In contrast, there was little liquid water content above 500 hPa from the MY scheme (Fig. 15e).

Figure 14.  Time-height cross-sections of domain-averaged (a) water vapor (q_v), (b) liquid (q_li = q_c + q_r), and (c) solid (q_so = q_i + q_s + q_g + q_h) hydrometeors from the simulation with the IMY scheme. The averaging area is along the rainband near the coastline, marked with a red rectangle in Fig. 10. The time-averaged profiles of the corresponding section are provided in panels (d) q_v, (e) q_li, and (f) q_so, respectively.

Figure 15.  Same as Fig. 14 but for the simulation with the MY scheme.

A considerable population of solid hydrometeors appears above 600 hPa levels with a peak value of over 0.45 g kg^–1. It is noteworthy that graupel was the dominant frozen hydrometeor owning to plentiful water vapor in the clouds. Under a rich water vapor condition, cloud ice particles grew rapidly, forming many graupel particles (McCumber et al., 1991; Krueger et al., 1995; Franklin et al., 2005). In comparison, low values of snow and hail content were found in the present case, which is consistent with Wen et al. (2006), which may be one of the fundamental features within clouds of severe rainfall over (sub)tropical regions. In summary, this simulation found considerable water vapor, rainwater, and graupel in the cloud compared to other hydrometeor categories. Solid hydrometeors appeared at upper levels above 600 hPa, and rainwater was mainly concentrated below 600 hPa. More specifically, graupel and rainwater were the major components of solid and liquid hydrometeors, respectively. Significant differences can be found in the distribution of solid hydrometeors between the IMY and MY schemes. Specifically, the MY scheme provided much larger solid hydrometeors than the IMY scheme. Combined with the surface rainfall amounts (Fig. 10) and the vertical profiles of the hydrometeors (Fig. 13), we may conclude that the MY scheme overproduced solid hydrometeors in the upper levels.

From the simulations, one can see that the first occurrence of liquid hydrometeors occurred nearly two hours before solid hydrometeors. This means that the precipitation started from warm rain microphysics processes at lower levels, and then ice phase processes were activated as a consequence. In warm clouds, a certain amount of time (usually termed gestation period) is required for cloud droplet coalescence to form large drops at the beginning of a rainfall event (Berry, 1968; Lin et al., 2005; Yin et al., 2015). This is consistent with the findings of previous studies (Lou et al., 2003; Lin et al., 2005) based on ideal simulations in which raindrops originally result from the autoconversion process of cloud droplets. With the development of the cloud, cold rain processes are activated. Ascending air, rich with water vapor and liquid droplets, is transported to upper levels where favorable conditions exist for ice particle growth. As ice particles fall into the region with above freezing temperatures, the melting of ice particles contributes considerably to the rainwater. As a result, severe rainfall formed as a result of both warm and cold microphysical processes (Yin et al., 2018). A large number of studies (e.g., Fovell and Ogura, 1988; McCumber et al., 1991; Tao et al., 2003; Gao et al., 2006) proposed that ice phase microphysical processes are more efficient to form large size drops than that of warm cloud processes. Consequently, overproduced solid hydrometeors were responsible for the overestimation of surface rainfall.

Based on the standard Milbrandt-Yau (MY) scheme in the Weather Research and Forecasting (WRF) model version 3.5, a modified double-moment bulk cloud microphysics scheme (named IMY) has been developed, aiming to improve cloud and precipitation simulation/forecasts over the East Asian monsoon region (EAMR). Numerous modifications have been developed according to long-term cloud-precipitation observations over the EAMR, as well as from other microphysics schemes. Most importantly, the shape parameters for raindrop and snow size distributions are diagnosed based on their respective shape-slope relationships. The shape parameter for cloud droplet size distribution depends on the total cloud number concentration. Though it is an additional computational skill to diagnose a shape parameter, it comes with little additional computational cost. A new approach for auto-conversion from cloud water to rainwater was also applied. In addition, the anomalously large ice nuclei concentration over East Asia is considered in the IMY scheme.

The IMY scheme has been incorporated into the WRF model and tested by modeling two heavy rainfall events (a mid-latitude extreme rainfall event and an extreme local rainfall event in the warm sector over the subtropical region) over the EAMR. Overall, the WRF model coupled with the IMY scheme worked very well in the triple-nested configuration with horizontal resolutions of 36 km, 12 km, and 4 km. In general, the WRF model with the IMY scheme can reproduce reliable results for the two selected heavy rainfall events. The simulations with the IMY scheme were comparable to those with the MY scheme and the observations. As for the mid-latitude extreme rainfall over northern China, the IMY scheme successfully reproduced the heavy rainfall with amounts over 200 mm in addition to their spatial distributions of cloud coverage and their evolution. Several embedded convective cells were successfully reproduced by the IMY scheme, while the organized convective cells were indiscernible in the simulation with the MY scheme. In the simulation of the local extreme rainfall in the warm sector over southern China, both the IMY and MY schemes can reproduce the major convective systems and the surface rainfall events in the warm sector over southern China, evidenced by comparing with the observations from rain gauges and ground-based weather radars. One of the major components of heavy rainfall along the coastal line from the southwest toward to northeast over southern China, associated with strong radar reflectivity in low levels of the convective system, was well reproduced by the IMY scheme. However, a spatial shift existed between simulation and observation. Generally speaking, considerable liquid water (cloud water) was allowed in the upper levels (above 500 hPa) in the IMY scheme, while there was little liquid water content in the MY scheme. In contrast, the MY scheme produced much larger solid hydrometeors than the IMY scheme. Considering the surface rainfall and vertical profiles of hydrometeors, we conclude that the overproduction of solid hydrometeors was responsible for the overestimation of rainfall by the MY scheme.

In addition to the simulations in this paper, the IMY scheme has been used to study a short squall line over the Yangtze River−Huaihe River (Li et al., 2014), a snow and rain mixing rainfall event over northern China (Zhang et al., 2016), and a heavy rainfall event over southern China (Yin et al., 2018). The preliminary evidence from the case studies indicates that the IMY scheme is potentially promising.

Although good results have been obtained from both the ideal and real-data case simulations, more tests are required to evaluate the performance of the IMY scheme because the precipitation over the EAMR is relatively complicated. For example, precipitation from the southwest vortex over the Sichuan Basin is likely the result of various processes associated with the Tibetan Plateau. Additionally, lasting rainfall over the Yangtze River is closely connected to the quasi-stationary Mei-Yu front. To date, the IMY scheme has only been tested by modeling several cases. In the near future, many experiments will be launched to further validate the IMY scheme for various weather events. Systematic experiments over long periods (at least a month) and detailed comparisons with multiple observations will be performed using the datasets from field experiments, such as the Southern China Monsoon Rainfall Experiment (SCMREX). Comparisons with other microphysics schemes are also needed. Moreover, to document microphysical processes, additional case studies, such as one including a budget analysis for each hydrometeor, should be accounted for in future research efforts (Huang et al., 2013).

Although good results have been achieved, it is still troublesome to explain the effects of cloud microphysical processes on the improvements because cloud microphysical processes and dynamical processes interact in complex ways (Gross et al., 2018). Consequently, considerable attention should be paid to interactions between cloud microphysical components and dynamical processes. Moreover, a Lagrangian particle-based approach may offer advantages for emerging techniques to model cloud microphysics (Grabowski et al., 2018).

Acknowledgements. This work is jointly supported by the National Natural Science Foundation of China (Grant No. 42075083), National Key Research and Development Program of China (Grant No. 2019YFC1510400), Guangdong Major Project of Basic and Applied Basic Research (Grant No. 2020B0301030004), and the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK010402). Particular thanks go to Jason Milbrandt (at the Meteorological Research Division/Environment and Climate Change, Canada) for using and modifying the code.