Mitigating hail overforecasting in the two-moment Milbrandt–Yau microphysics scheme (v2.25.2_beta_04) in WRF (v4.5.1) by incorporating the graupel spongy wet growth process (MY2_GSWG v1.0)

Hail forecasting using numerical models remains a challenge due to the uncertainties and deficiencies in microphysics schemes. In this study, we assessed the hail simulation performance of the two-moment Milbrandt–Yau (MY2) microphysics scheme within the Weather Research and Forecast (WRF) model by simulating three heavy rainfall events in Meiyu systems in which hail was rarely observed, rather than focusing on hail cases as done in previous research. Simulation results showed that the MY2 scheme produced noticeable hail in these rainstorms. Further analysis revealed that the overprediction of hail was caused by the imperfect graupel-to-hail conversion parameterization method adopted in the MY2 scheme. By incorporating the graupel spongy wet growth process, the modified scheme significantly mitigated the hail overforecasting. Moreover, the modified MY2 scheme kept the ability to simulate hail in real hail cases, demonstrating its ability to differentiate between heavy rainfall and hail events – a distinction the original scheme lacked. By comparing simulations of both rainstorms and hailstorms, it is concluded that the upward transport of large raindrops near the 0 °C level is critical for the graupel-to-hail conversion.

1 Introduction

Hailstorms are high-impact weather events with significant damage potential, amounting to billions of US dollars in losses annually around the world (Zhang et al., 2008; Changnon, 2009; Brown et al., 2015; Kunz et al., 2018). However, accurately predicting hail remains a challenge in current numerical weather prediction models due to the need for precise representations of environmental conditions, the dynamic structure of hailstorms, and the associated microphysical processes (Wellmann et al., 2020).

Traditional hail forecasting methods primarily rely on environmental parameters derived from soundings to predict the likelihood of severe thunderstorms, rather than explicitly forecasting the conditions favorable for hail formation and development (Allen et al., 2011). Moreover, predicting hail size using sounding-derived parameters is particularly challenging, with limited success in previous studies (Smith et al., 2012; Tuovinen et al., 2015). Hailstones are influenced not only by large-scale environmental conditions but also by the intensity of the parent hailstorm, the structure of the wind field, and the microphysical processes involved in hail development. These processes typically occur on small spatial and temporal scales, while sounding data represent large-scale, time-averaged meteorological conditions. As a result, these methods often fail to accurately capture the complex mechanisms underlying hail formation and development (Allen et al., 2019).

In recent years, the application of numerical models in hail forecasting has gained increasing attention. The 1D hail growth scheme, called the HAILCAST scheme (Adams-Selin and Ziegler 2016), has shown reasonable accuracy in predicting hail size (Gagne et al., 2017; Adams-Selin et al., 2019, 2023). However, this scheme does not directly simulate the formation of hail or the growth of hail populations. Instead, it places five hail embryos at predefined vertical levels in suitable locations once some specific dynamic conditions are met and then assesses how large these hailstones will grow through a series of microphysical processes. Since this scheme does not calculate microphysical processes associated with hail populations, it does not provide feedback on the dynamic, thermal, and other hydrometeor fields. Furthermore, accurate simulation of hail diameter in different regions and cases depends on parameter tuning within the HAILCAST scheme (Adams-Selin et al., 2019).

Using microphysical schemes that can simulate hail populations (referred to as the hail microphysical scheme in this study, or HMS) is a more reasonable and comprehensive approach to simulating hail cases, especially in terms of the mutual feedback between thermodynamics and cloud microphysics. Previous studies have demonstrated that incorporating hail and associated microphysical processes could improve the simulation of storm structure, cold pool intensity, and radar echo morphology (Bryan and Morrison, 2012; Tao et al., 2016). Employing HMSs allows for forecasting hail size, investigating the mechanisms of hail formation, and exploring the impact of environmental conditions on hail formation and evaluation (e.g., Milbrandt and Yau, 2006a; Dennis and Kumjian, 2017; Li et al., 2017; Yin et al., 2021; Wu et al., 2022; Dong et al., 2024).

Before these microphysics schemes can be reliably used, their hail simulation performance must be evaluated against observations. Previous research has often evaluated the accuracy of hail simulation by simulating specific hail cases and comparing model-diagnosed hail size and hail coverage with those derived from operational radar data through the maximum estimated size of hail (MESH) algorithm (Witt et al., 1998) or hydrometeor classification algorithms (HCAs; Park et al., 2009; Ortega et al., 2016). For example, Luo et al. (2017, 2018) simulated two hail events in eastern China using different configurations of the Milbrandt scheme (Milbrandt and Yau, 2005a, b) and found that the three-moment Milbrandt scheme yielded surface hail diameters and hail coverage areas that more closely matched observations. However, they also noted that the Milbrandt scheme inaccurately simulated larger hail diameters below the freezing level. Labriola et al. (2019a) demonstrated that this overestimation of hail diameter is due to the continued growth of hail by collecting cloud droplets and raindrops below the melting layer. Additionally, Labriola et al. (2017, 2019a, b) compared observed and simulated hail diameters across multiple hail events and identified various issues with multiple HMSs. For example, the National Severe Storms Laboratory (NSSL) scheme (Mansell et al., 2010) failed to simulate large hail, the Milbrandt scheme showed excessive melting at lower levels, and the P3 scheme (Morrison and Milbrandt, 2015) simulated high-density ice particles only near the surface, while low-density ice particles dominated above the melting layer.

However, there are some issues with quantitatively evaluating the model's ability to simulate hail size based on the observations. First, methods like MESH and HCA just estimate hail diameter based on radar data, which are not direct measurements of hail size. Second, the HMS does not directly forecast hail diameter; instead, it diagnoses the maximum hail diameter using a method that derives it from the predicted mass and number mixing ratios of hail particle populations. The parameter in this diagnostic method is artificially set, introducing uncertainties in the diagnosed hail diameter. Furthermore, numerous microphysical processes can influence the simulated hail diameter, making it challenging to pinpoint which specific process is responsible when the simulation produces an inaccurate hail diameter.

A more effective and fundamental approach to evaluate HMSs is to qualitatively assess the presence or absence of hail. Binary classification of hail presence based on HCA is more reliable than estimated hail size. Moreover, focusing on the presence or absence of hail when evaluating microphysics schemes shifts the emphasis to processes associated with hail formation, making it easier to identify the sources of simulation biases. We believe that a robust HMS should at least reliably distinguish between the presence and absence of hail, accurately producing hail in cases where it forms and not producing it in cases where it does not.

In this study, we evaluate the capability of a selected HMS (MY2 scheme) to simulate the presence or absence of hail using data from heavy rainfall events within the Meiyu system in eastern China. Previous research has shown that in typical Meiyu systems characterized by high total precipitable water and uniformly moist conditions, heavy rainfall is more likely to occur, while hail is rarely observed (Chen et al., 2019; Wang et al., 2019; Yao et al., 2020). These rainstorms with no hail provide excellent examples for examining the ability of HMS to simulate hail. The description of the model configuration, cases, and data used in this study is provided in the next section. Section 3 details the model evaluation results, modifications to the selected hail microphysics scheme, and reevaluation of the modified scheme using data from both heavy rainfall and hail events. Section 3 also includes an interpretation of why hail forms in hail events but does not form in heavy rainfall events. Finally, conclusions and a discussion regarding the modified HMS are presented in Sects. 4 and 5, respectively.

2 Model and methods

2.1 Model configuration

This study employs the Advanced Weather Research and Forecasting (WRF; Skamarock et al., 2021) model version 4.5.1 to simulate three heavy rainfall events and three hail events that occurred in eastern China. The 0.25° × 0.25° 1 h interval ERA5 reanalysis data are used as the initial and boundary conditions for each case. The digital filter initialization method is adopted to remove initial model imbalance. Simulations are conducted with three two-way nested domains, with the horizontal grid spacings being 9, 3, and 1 km, respectively. The time step for the innermost domain is 2 s. All the simulation analyses are based on the results provided by the innermost 1 km spacing domain. Details about the model configuration for each case are shown in Table 1. The Noah Multiphysics land surface model (Niu et al., 2011) and the Rapid Radiative Transfer Model for GCMs (RRTMG) radiation scheme (Pincus et al., 2003) are used for all simulations. The Mellor–Yamada–Janjic (MYJ) planetary boundary layer scheme (Janjić, 1994) is used when simulating heavy rainfall events, whereas the Asymmetric Convective Model version 2 (ACM2) scheme (Pleim, 2007) was chosen for hail simulations due to its more accurate representation of the timing and location of convective initiation compared to observations.

Table 1 Model configurations.

This study selected the MY2 scheme for evaluating hail simulations due to its recognition as a classic bulk microphysics scheme known for its ability to forecast hail. This scheme has significantly influenced the development of other microphysics schemes. For example, the graupel-to-hail conversion parameterization from the MY2 scheme, as mentioned in this study, has been adopted and further refined by the NSSL, National Taiwan University (NTU, Tsai and Chen, 2020), and spectral bin microphysics (SBM, Khain et al., 2010a) schemes. Following Labriola et al. (2019a), the accretion of hail by collecting liquid water beneath the melting layer is prohibited in this study to avoid the unreasonable increase in hail size within this layer.

In addition to surface hail reports, the observed hail occurrence facts are also obtained from the hydrometeor classification algorithm (HCA) applying to S-band polarimetric radar. The HCA in this study comes from Dolan et al. (2013), which is a fuzzy logic algorithm based on the construction of membership functions for each hydrometeor type using multiple polarimetric radar variables, including radar reflectivity (Z_H), differential reflectivity (Z_DR), specific differential phase (K_DP), and copolar correlation coefficient (ρ_hv). The dominant hydrometeor type in a given radar sample volume is identified (Vivekanandan et al., 1999). Although 10 types of hydrometeors are identified by the HCA, whether the identified type is hail is the primary concern of this paper. Hail and graupel are primarily distinguished by the combined distribution ranges of multiple polarimetric radar variables. For example, the radar reflectivity factor (Z_H) for graupel is typically in the range of 34.1–54.5 dBZ, while for hail, it is between 48–76.6 dBZ. The differential reflectivity (Z_DR) for graupel is usually between 0.4–2.8 dB, whereas for hail, it is between −0.4 and 0.7 dB.

2.2 Modification of the MY2 scheme

According to the MY2 scheme, graupel particles that have just converted to hail through wet growth are envisioned as ice balls with a loose ice core with a bulk density of 400 kg m⁻³, covered by a thin film of unfrozen water. However, real hailstone images show that even small hailstones, which primarily undergo dry growth and lack transparent layers, have a hard and high-density ice core. Graupel particles have internal voids, resulting in a bulk density lower than that of pure ice. During wet growth, unfrozen supercooled water first infiltrates these voids and freezes, increasing the bulk density of the graupel. Then, the remaining unfrozen water forms a film on the surface. This process is known as spongy wet growth (Pruppacher and Klett, 1997; Khain and Pinsky, 2018). Most microphysics schemes do not account for this process. Filling the internal voids requires more supercooled water, which can prevent graupel particles from converting to hail too easily near the 0 °C level; this will be shown in the following analysis.

Furthermore, larger particles collect more supercooled water upon collision, making it more likely to undergo wet growth. Therefore, when the amount of supercooled water is insufficient to activate wet growth for all the graupel particles with different sizes, only those larger than a certain size will undergo wet growth, while smaller particles continue to undergo dry growth. In the modified MY2 scheme, graupel particles are allowed to convert to hail only under specific conditions. This conversion occurs when the total quantity of supercooled cloud droplets and raindrops collected by particles larger than a critical size exceeds the quantity required for their wet growth and for filling their internal voids to reach a density of 900 kg m⁻³ within one integration time step. We use the incomplete gamma function to describe the integration over the graupel spectrum from a certain particle size D^* to infinity. The total quantity of supercooled cloud droplets collected by graupel particles larger than D^* per second Q_cg (unit: kg m⁻³ s⁻¹) is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{Q}_{\mathrm{cg}}& ={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\underset{D^{\ast }}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left|V_{\mathrm{g}}\left(D_{\mathrm{g}}\right)-V_{\mathrm{c}}\left(D_{\mathrm{c}}\right)\right|{\left(D_{\mathrm{g}}+D_{\mathrm{c}}\right)}^{\mathrm{2}}\\ & m_{\mathrm{c}}\left(D_{\mathrm{c}}\right)E_{\mathrm{cg}}{N}_{\mathrm{g}}\left(D_{\mathrm{g}}\right)N_{\mathrm{c}}\left(D_{\mathrm{c}}\right)\mathrm{d}{D}_{\mathrm{c}}\mathrm{d}{D}_{\mathrm{g}},\end{array}\end{array}$$

where D_g (D_c), V_g (V_c=0), and N_g (N_c) are the diameter, terminal fall velocity, and size distribution of graupel (cloud), respectively. $V_{\mathrm{g}}=(\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{air}}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}{a}_{\mathrm{g}}{D}_{\mathrm{g}}^{b_{\mathrm{g}}}$, where ρ_air is air density and ρ₀=1.225 kg m⁻³. V_c=0. $m_{\mathrm{c}}=\frac{\mathit{\pi }}{\mathrm{6}}{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{c}}^{\mathrm{3}}$ is the mass of a cloud droplet with diameter D_c, where ρ_w is water density. E_cg is the collecting efficiency between graupel and cloud, which is a function of the mass-weighted mean diameters of cloud droplets and graupel particles (Cober and List, 1993). In the Milbrandt scheme, cloud size distribution is parameterized with $N_{\mathrm{c}}\left(D\right)=N_{\mathrm{0}\mathrm{c}}{D}^{{\mathit{\alpha }}_{\mathrm{c}}}{e}^{-{\mathit{\lambda }}_{\mathrm{c}}D}$, where $N_{\mathrm{0}\mathrm{c}}=N_{\mathrm{Tc}}\frac{\mathrm{1}}{\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}})}{\mathit{\lambda }}_{\mathrm{c}}^{\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}}}$, where N_Tc is the total number concentration and α_c=1. By substituting these equations and integrating them over the cloud and graupel spectrum, we can obtain the following:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{Q}_{\mathrm{cg}}& ={\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{air}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}{a}_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{24}}}{E}_{\mathrm{cg}}{\displaystyle \frac{N_{\mathrm{Tc}}{N}_{\mathrm{0}g}}{\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}})}}\\ & \left[{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{3}+b_{\mathrm{g}},{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast })\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}}+\frac{\mathrm{3}}{v_{\mathrm{c}}})}{{\mathit{\lambda }}_{\mathrm{c}}^{\mathrm{3}}}}\right.\\ & +{\displaystyle \frac{\mathrm{2}\mathrm{\Gamma }(\mathrm{2}+b_{\mathrm{g}},{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast })\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}}+\frac{\mathrm{4}}{v_{\mathrm{c}}})}{{\mathit{\lambda }}_{\mathrm{c}}^{\mathrm{4}}}}\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{1}+b_{\mathrm{g}},{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast })\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_{\mathrm{c}}+\frac{\mathrm{5}}{v_{\mathrm{c}}})}{{\mathit{\lambda }}_{\mathrm{c}}^{\mathrm{5}}}}\right].\end{array}\end{array}$$

Similarly, the total quantity of supercooled raindrops collected by graupel particles larger than D^* per second Q_rg (unit: kg m⁻³ s⁻¹) is

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{Q}_{\mathrm{rg}}& ={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\underset{D^{\ast }}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left|V_{\mathrm{g}}\left(D_{\mathrm{g}}\right)-V_{\mathrm{r}}\left(D_{\mathrm{r}}\right)\right|{\left(D_{\mathrm{g}}+D_{\mathrm{r}}\right)}^{\mathrm{2}}\\ & m_{\mathrm{r}}\left(D_{\mathrm{r}}\right)E_{\mathrm{rg}}{N}_{\mathrm{g}}\left(D_{\mathrm{g}}\right)N_{\mathrm{r}}\left(D_{\mathrm{r}}\right)\mathrm{d}{D}_{\mathrm{r}}\mathrm{d}{D}_{\mathrm{g}},\end{array}\end{array}$$

according to Murakami (1990). $\left|v_{\mathrm{g}}\left(D_{\mathrm{g}}\right)-v_{\mathrm{r}}\left(D_{\mathrm{r}}\right)\right|\approx \sqrt{(V_{Q\mathrm{g}}-V_{Q\mathrm{r}}{)}^{\mathrm{2}}+\mathrm{0.04}{V}_{Q\mathrm{g}}{V}_{Q\mathrm{r}}}$, where V_Qg and V_Qr are the mass-weighted mean fall velocity of graupel and rain, respectively. Rain size distribution is also parameterized as $N_{\mathrm{r}}\left(D\right)=N_{\mathrm{0}\mathrm{r}}{D}^{{\mathit{\alpha }}_{\mathrm{r}}}{e}^{-{\mathit{\lambda }}_{\mathrm{r}}D}$, where α_r=0. Thus, we can get

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{Q}_{\mathrm{rg}}& =\sqrt{(V_{Q\mathrm{g}}-V_{Q\mathrm{r}}{)}^{\mathrm{2}}+\mathrm{0.04}{V}_{Q\mathrm{g}}{V}_{Q\mathrm{r}}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{24}}}\\ & E_{\mathrm{rg}}{N}_{\mathrm{Tr}}{N}_{\mathrm{0}g}{e}^{-{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast }}\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{4}\right)}{{\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{3}}}}\left({\displaystyle \frac{D^{\ast \mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{2}{D}^{\ast }}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{3}}}}\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{2}\mathrm{\Gamma }\left(\mathrm{5}\right)}{{\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{4}}}}\left({\displaystyle \frac{D^{\ast }}{{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}\right)+{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{6}\right)}{{\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{5}}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{g}}}}\right].\end{array}\end{array}$$

Based on Musil (1970), Lin et al. (1983), and Milbrandt and Yau (2005a), the critical amount of supercooled water required for graupel particles larger than D^* to undergo wet growth Q_gwet (unit: kg m⁻³ s⁻¹) is

$$\begin{array}{}\text{(5)}& Q_{\mathrm{gwet}}=\underset{D^{\ast }}{\overset{\mathrm{\infty }}{\int }}\mathrm{2}\mathit{\pi }{D}_{\mathrm{g}}\mathrm{Vent}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{L}_{\mathrm{v}}{D}_{\mathrm{v}}\mathrm{\Delta }\mathit{\rho }-K_{\mathrm{a}}{T}_{\mathrm{c}}}{L_{\mathrm{f}}+C_{\mathrm{w}}{T}_{\mathrm{c}}}}{N}_{\mathrm{g}}\left(D_{\mathrm{g}}\right)\mathrm{d}{D}_{\mathrm{g}},\end{array}$$

where L_v is the latent heat of vaporization, D_v is the molecular diffusion coefficient of water, K_a the thermal conductivity of air, T_c the Celsius temperature, and C_w the specific heat of water. Vent is the bulk ventilation coefficient for graupel, which is parameterized as $\mathrm{Vent}=\mathrm{0.78}+\mathrm{0.308}S{c}^{\frac{\mathrm{1}}{\mathrm{3}}}R{e}^{\frac{\mathrm{1}}{\mathrm{3}}}$, where Sc is the Schmidt number and Re the Reynolds number. $Re=\frac{D_{\mathrm{g}}{V}_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{air}}}{\mathit{\mu }}$, where μ is the dynamic viscosity of air. In the derivations of Eq. (5) by Musil (1970) and Milbrandt and Yau (2005a), graupel particles were assumed to be spherical, and their capacitance was set equal to the graupel particle radius. Equation (5) indicates that when the supercooled water content exceeds Q_gwet, all graupel particles with diameters larger than D^* will undergo wet growth. By substituting these equations and integrating them over the graupel spectrum from particles larger than D^*, we can obtain the following:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{Q}_{\mathrm{gwet}}& =\mathrm{2}\mathit{\pi }{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{L}_{\mathrm{v}}{D}_{\mathrm{v}}\mathrm{\Delta }\mathit{\rho }-K_{\mathrm{a}}{T}_{\mathrm{c}}}{L_{\mathrm{f}}+C_{\mathrm{w}}{T}_{\mathrm{c}}}}\\ & N_{\mathrm{0}\mathrm{g}}\left[{\displaystyle \frac{\mathrm{0.78}{e}^{-{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast }}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+\mathrm{0.308}S{c}^{\frac{\mathrm{1}}{\mathrm{3}}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}}{\mathit{\mu }}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{air}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}{a}_{\mathrm{g}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right.\\ & \left.{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2.5}+\frac{b_{\mathrm{g}}}{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast }\right)}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2.5}+\frac{b_{\mathrm{g}}}{\mathrm{2}}}}}\right].\end{array}\end{array}$$

The supercooled water mass required to fill the internal voids of graupel larger than D^* per time step Q_fill (unit: kg m⁻³) is

$$\begin{array}{}\text{(7)}& Q_{\mathrm{fill}}=\underset{D^{\ast }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}({\mathit{\rho }}_{\mathrm{h}}-{\mathit{\rho }}_{\mathrm{g}})D_{\mathrm{g}}^{\mathrm{3}}{N}_{\mathrm{g}}\left(D_{\mathrm{g}}\right)\mathrm{d}{D}_{\mathrm{g}},\end{array}$$

where ρ_h and ρ_g are the bulk densities of graupel and hail, respectively. In the MY2 scheme, they are predefined as ρ_h=900 kg m⁻³ and ρ_g=400 kg m⁻³. Integrating Eq. (7) we can get

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{Q}_{\mathrm{fill}}& ={\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}({\mathit{\rho }}_{\mathrm{h}}-{\mathit{\rho }}_{\mathrm{g}})N_{\mathrm{0}\mathrm{g}}{e}^{-{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast }}\\ & \left({\displaystyle \frac{D^{\ast \mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{3}{D}^{\ast \mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{6}{D}^{\ast }}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{6}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{4}}}}\right).\end{array}\end{array}$$

In the modified MY2 scheme, we specify that graupel particles larger than D^* undergo wet growth and convert to hail only when the amount of supercooled liquid water they collected ((Q_rg+Q_cg)dt) exceeds the sum of the amount required to raise their surface temperature to 0 °C (Q_gwetdt) and the amount needed to fill their internal voids (Q_fill). That is,

$$\begin{array}{}\text{(9)}& \left(Q_{\mathrm{rg}}+Q_{\mathrm{cg}}\right)\mathrm{d}t\ge Q_{\mathrm{gwet}}\mathrm{d}t+Q_{\mathrm{fill}},\end{array}$$

where dt is the model integration time step. D^* ranges from 0.5 to 25 mm, with a 0.5 mm interval, and the inequality condition in Eq. (9) is checked for each value. When D^* satisfies the condition but (D^* – 0.5 mm) does not, the final D^* is determined through linear interpolation. The mass and concentration of particles larger than D^* in the graupel spectrum are then transferred to the hail category. The mass of graupel that converted to hail is calculated as that in the original MY2 scheme, which is

$$\begin{array}{}\text{(10)}& Q_{\mathrm{cngh}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{\mathit{\rho }}_{\mathrm{g}}{N}_{\mathrm{0}g}{e}^{-{\mathit{\lambda }}_{\mathrm{g}}{D}^{\ast }}\left({\displaystyle \frac{D^{\ast \mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{3}{D}^{\ast \mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{6}{D}^{\ast }}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{6}}{{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{4}}}}\right).\end{array}$$

The modified Milbrandt scheme is referred to as MY2_Mod, while the original one is referred to as MY2_Ori.

2.3 Case overview

We selected three typical heavy rainfall events within a Meiyu front system and three hail events that occurred in eastern China between 2021 and 2024 to evaluate the hail simulation performance of the MY2 scheme. The heavy rainfall events took place on 23–24 June 2022 (Rainstorm1), 17–18 June 2023 (Rainstorm2), and 6–7 July 2023 (Rainstorm3). As shown in Fig. 1, the synoptic conditions of Rainstorm1 and Rainstorm3 are quite similar: they are both located in relatively flat westerly belts, situated south of the 500 hPa trough and to the northwest of the western Pacific subtropical high. These two heavy rainfall events are characterized by abundant water vapor, with the precipitable water exceeding 64 and 56 mm, respectively. Combined with low-level southwest moisture transport, the 12 h cumulative precipitation reached a peak of 189.9 mm in Rainstorm1 (from 12:00 on 23 June to 00:00 on 24 June 2022, UTC) and 235.1 mm in Rainstorm3 (from 18:00 on 6 July to 06:00 on 7 July 2023), with the maximum hourly rainfall reaching 118.8 and 87.3 mm h⁻¹, respectively. Rainstorm2 also occurred in a relatively flat 500 hPa westerly belt to the northwest of the western Pacific subtropical high, but the synoptic-scale forcing to its north was less pronounced than in the other two cases, and the precipitable water was also slightly lower, mostly ranging between 48–56 mm. The maximum hourly rainfall in this event reached 64.6 mm, with 12 h cumulative precipitation of 149.8 mm (from 18:00 on 17 July to 06:00 on 18 July 2023).

The three selected hail events occurred on 30 April 2021 (Hailstorm1), 24 June 2022 (Hailstorm2), and 25 March 2024 (Hailstorm3). The synoptic conditions of Hailstorm1 and Hailstorm2 are similar, with both located in the northwesterly flow behind the 500 hPa trough and characterized by relatively low precipitable water, which was below 16 and 24 mm, respectively. In Hailstorm1, a supercell formed in southern Shandong Province at 07:00 UTC (15:00 LST) on 30 April 2021 and gradually developed into a squall line, moving from northwest to southeast. This hail event brought strong winds (with maximum gusts reaching 45.4 m s⁻¹), hail (with a maximum diameter of 3–5 cm), and moderate rainfall (with a maximum hourly rainfall of 37.5 mm) to Shandong and Jiangsu provinces. In Hailstorm2, several supercells formed in southern Shandong Province at 06:00 UTC (14:00 LST) on 24 June 2022 and then gradually moved southward, bringing short-term heavy rainfall (with a maximum hourly rainfall of 61.5 mm) and large hailstones with diameters of 5–8 cm. In contrast to these two hail events, which occurred in the northwesterly flow behind the 500 hPa trough, Hailstorm3 developed in the southwesterly flow ahead of a short trough at 500 hPa, with more precipitable water (40–48 mm). In this event, a supercell formed in western Zhejiang Province at 06:00 UTC (14:00 LST) and gradually developed into a squall line, affecting the entire province from west to east. This hail event produced short-term heavy rainfall (with a maximum hourly rainfall of 49.9 mm) and large hailstones with maximum diameters of 5–6 cm.

Figure 1 Synoptic conditions obtained from the ERA5 reanalysis at 12:00 UTC on (a) 23 June 2022 (Rainstorm1), (b) 17 June 2023 (Rainstorm2), and (c) 6 July 2023 (Rainstorm3) and at 00:00 UTC on (d) 30 April 2021 (Hailstrorm1), (e) 24 June 2022 (Hailstrorm2), and (f) 24 March 2024 (Hailstrorm3). Color shading represents precipitable water (unit: mm). Black solid lines represent geopotential height (unit: gpm). White arrows represent the wind field at 850 hPa. Pink boxes in each panel represent the location of the two innermost domains.

3 Results

3.1 Results with the original MY2 scheme

Figure 2 presents a comparison between observed and simulated column maximum radar reflectivity for three heavy rainfall events. The observed radar echoes for these events exhibit an east–west orientation, presenting a typical feature of a Meiyu front system. The model successfully reproduces this east–west distribution of radar echoes. In addition, the location of the simulated radar echoes closely aligns with the observations. This indicates that the numerical model successfully captures the main characteristic of these three Meiyu front precipitation events, making the simulation results reliable for further analysis.

Figure 2 Comparison of the observed (a, b, c) and simulated (d, e, f) column maximum radar reflectivity at (a, d) 17:20 UTC on 23 June 2022, (b, e) 20:40 UTC on 17 June 2023, and (c, f) 18:00 UTC on 6 July 2023. The “×” symbol indicates the presence of at least one point in the vertical direction where the reflectivity factor of hail exceeds the sum of all other hydrometeors. Black lines represent the locations of the vertical cross-sections shown in Fig. 3. Gray lines represent provincial boundaries.

As the observation radar used is a dual-polarization radar, its polarimetric radar variables, such as radar reflectivity, differential reflectivity (Z_DR), and specific differential propagation phase (K_DP), could be used for hydrometeor classification. The first row in Fig. 2 shows the 2D classification result, with “×” indicating the presence of identified hailstones in the vertical direction at that position. However, as shown in the figure, no hail was identified in all three heavy rainfall events.

It should also be noted that the absence of hail, as mentioned above, refers to the non-detection of hail using dual-polarization radar data. This does not imply that there was not a single hail particle or hail embryo in the real convective clouds. It just indicates that, in these three events, the convective clouds did not contain enough large or numerous hailstones that could be detected by dual-polarization radar. In the MY2 scheme, the hail category represents ice-phase particles with higher density and faster terminal velocities, which are not entirely equivalent to the hail identified through HCA in observations. To align with the hail identification criteria used in observations, we stipulate that simulated hail in the model is considered to exist at a model grid only when it meets the following three conditions: (1) its mass mixing ratio exceeds 0.1 g kg⁻¹, (2) its mass-weighted mean diameter exceeds 2 mm and is larger than the graupel mass-weighted mean diameter at the same model grid, and (3) its equivalent reflectivity factor (Z_e) must exceed the combined Z_e of all other hydrometeor species (cloud, rain, ice, snow, and graupel). It means that only points where hailstones have a non-negligible mass, are larger than graupel particles, and contribute the most to the total radar reflectivity are identified as positions where hail is present. This method is in principle similar to, though not identical to, the hail identification algorithm adopted in observations. We did not use the method of first calculating the dual-polarization variables of the simulation results with a radar simulator and then identifying hail using the same hydrometeor classification algorithm as in the observations, since the MY2 scheme lacks proper treatment of ice particle shapes, which are essential for polarization sensitivity, resulting in limited benefits from using a simulator (Matsui et al., 2020).

The second row in Fig. 2 illustrates that the simulations using the original MY2 scheme significantly overestimated hail amounts. Hail was present in most areas where the column maximum reflectivity exceeded 40 dBZ, contrary to observations that detected no hail. Figure 3 illustrates the ratio of hail occurrence in regions where the column maximum reflectivity exceeds 40 dBZ, comparing observations with simulations across different periods of three heavy rainfall events. In the model, hail was detected in approximately 40 % of these regions, with this ratio remaining consistent for 12 h during all three events. In contrast, observations showed that this ratio remained near zero for extended periods. The persistent overestimation of hail across multiple cases and over prolonged durations suggests that the MY2 scheme's tendency to overpredict hail is a systematic and recurring issue, rather than an isolated occurrence in specific moments or cases.

Figure 3 Comparison of hail percentages in areas where column maximum radar reflectivity exceeds 40 dBZ between observations (solid lines) and simulations using the original (dashed lines) and modified (dotted lines) Milbrandt two-moment schemes. Red, blue, and orange lines correspond to Rainstorms 1, 2, and 3, respectively. The red, blue, and orange x-axis labels represent time spans from 12:00 UTC on 23 June to 00:00 UTC on 24 June, 16:00 UTC on 17 June to 04:00 UTC on 18 June, and 15:00 UTC on 6 July to 03:00 UTC on 7 July.

Figure 4 further compares the observed and simulated radar reflectivity and the spatial distribution of graupel and hail particles in the vertical direction. Similar to the previous method, regions identified as graupel in the observations are marked with dots. In the simulations, dots represent regions where the graupel mass mixing ratio exceeds 0.1 g kg⁻¹ and its equivalent reflectivity factor exceeds the sum of all the other hydrometeors (cloud, rain, ice, snow, and hail). Although the simulated radar echoes resemble the observed ones structurally, some discrepancies are present. Observed radar reflectivity values exceeding 40 dBZ are mainly concentrated below 6 km. In contrast, in the simulations, while most regions with reflectivity over 40 dBZ are also below 6 km, echoes exceeding 40 dBZ are occasionally found above this level. Significant differences are also noted in the vertical distribution of graupel and hail. In the observations, no hail is detected in the convection, and graupel particles are primarily located above 5000 m (near the 0 °C level) in regions with reflectivity of 30 dBZ or higher. Even in deep moist convection, where the 30 dBZ echo top height reached nearly 12 km (Fig. 3c), no hail was detected. In contrast, the simulations show that while graupel is present in the upper parts of the convection, a substantial amount of hail appears above the 0 °C level, forming several column-shaped regions with reflectivity exceeding 40 dBZ.

Figure 4 Comparison of the cross-sections of (a, b, c) observed and (d, e, f) simulated radar reflectivity at (a, d) 17:20 UTC on 23 June 2022, (b, e) 20:40 UTC on 17 June 2023, and (c, f) 18:00 UTC on 6 July 2023. The “×” symbols (dots) indicate that the reflectivity factor of hail (graupel) exceeds the sum of all other hydrometeors at that location. Gray lines represent the 0 °C level.

Figure 5 further illustrates the reflectivity, mass mixing ratio, mass-weighted mean diameter, and particle size distribution of graupel and hail in three heavy rainfall cases. The reflectivity of both graupel and hail is calculated according to the MY2 scheme based on the predicted mass and concentration mixing ratios. In these three heavy rainfall cases, graupel predominantly occurs above the 0 °C layer, with mass mixing ratios reaching up to 7 g kg⁻¹. In regions where the its mass mixing ratio exceeds 1 g kg⁻¹, the corresponding reflectivity typically ranges from 25 to 45 dBZ. The mass mixing ratio of hail is smaller than that of graupel, peaking at only 3 g kg⁻¹ in the figure. However, its radar reflectivity often exceeds 40 or 45 dBZ. Hail with a mass mixing ratio greater than 0.1 g kg⁻¹ is observed only above 3000 m, indicating that hail is rare in the lower layers and at the surface. Based on the criteria set in this study, it can be concluded that hail does not reach the surface in the model, despite the fact that the reflectivity of hail consistently extends from aloft to the surface, with values exceeding 40 dBZ. The hail reflectivity below 4 km may arise from the excessive size sorting of hail particles, as suggested by Milbrandt and Yau (2006b) and Milbrandt et al. (2021). Notably, in some regions of the three cases, hail unexpectedly appears below the 0 °C layer, while there is no hail present above the 0 °C layer, where only graupel exists. This suggests that in the model, hail may form from graupel as it passes through the 0 °C layer, which contradicts the widely accepted hail formation mechanisms.

The third column presents the mass-weighted mean diameter of graupel (D_mg) and hail (D_mh). In regions identified as hail in Fig. 4, D_mh is generally above 3 mm, while graupel has a mass-weighted diameter of around 2 mm. It indicates that hailstones are larger than graupel particles. The presence of larger hailstones is also evident in the particle size distributions (PSDs) shown in the fourth column. Although the total concentration of hail particles is much lower than that of graupel, the spectrum width of hail is larger. The concentration of hail particles larger than 5 mm at the three selected points is 4.86, 1.76, and 1.09 m⁻³, while the concentration of graupel particles larger than 5 mm does not exceed 0.2 m⁻³. Figure 5 indicates that although hail does not reach the surface in the three heavy rainfall cases in the model, there is indeed a sufficiently large number of hail particles with diameters exceeding those of graupel in the regions aloft that are identified as containing hail, further supporting the notion that the overforecasting of hail in the MY2 scheme is a real issue.

Figure 5 Recalculated radar reflectivity (color shading, units: dBZ) of graupel (a, e, i) and hail (b, f, j), mass-weighted mean diameter (c, g, k, units: mm) of graupel (gray solid line) and hail (color shading), and particle size distribution (d, k, l) of graupel (gray line) and hail (red line) for Rainstorms 1 (top row), 2 (middle row), and 3 (bottom row). Red stars in the third column indicate the locations of the plotted PSDs in the fourth column. Gray dashed lines in the first three columns represent the 0 °C layer.

The overprediction of hail also leads to a misunderstanding of the mechanisms of precipitation formation in the Meiyu system. Figure 6 presents the contribution of various microphysical processes to rainwater in regions where the rain rate exceeds 20 mm h⁻¹ for the three heavy rainfall cases. The primary source of rainwater in all three cases is the accretion of rainwater by collecting cloud droplets, with median contributions of 51.1 % (Rainstorm1), 60.6 % (Rainstorm2), and 62.8 % (Rainstorm3). This is reasonable, as these cases all have relatively high precipitable water (as shown in Fig. 1), and the warm layer with temperature higher than 0 °C is sufficiently deep (over 5 km). However, in all three cases, the median contribution of hail melting to rainwater is 22.5 % (Rainstorm1), 17.7 % (Rainstorm2), and 17.5 % (Rainstorm3), which is comparable to the contribution from graupel melting (20.6 %, 18.3 %, and 17.8 %). This suggests that, in the model, hail melting is the second-largest contributor to rainwater. This is in stark contrast to the observations, where hail was not detected at all.

Figure 6 Contribution of various microphysical processes to rainwater in regions where the rain rate exceeds 20 mm h⁻¹ for the three heavy rainfall cases. “Auto” stands for the autoconversion of cloud droplets to rainwater, “accr” represents the accretion of rain by collecting cloud droplets, and “sml”, “gml”, and “hml” are the melting of snow, graupel, and hail to form rain, respectively. The red, blue, and orange x-axis labels represent time spans from 12:00 UTC on 23 June to 00:00 UTC on 24 June, 16:00 UTC on 17 June to 04:00 UTC on 18 June, and 15:00 UTC on 6 July to 03:00 UTC on 7 July.

Where do these overpredicted hailstones originate? We investigated the hail formation processes in the MY2 scheme, including the graupel-to-hail conversion process, raindrop freezing process, and three-component freezing of raindrops that lead to hail formation. The result shows that, in terms of mass mixing ratio, the conversion of graupel to hail is the dominant pathway for hail formation. In the three heavy rainfall cases, this process accounted for 99.8 %, 98.1 %, and 99.9 % of hail formation (not shown). In the MY2 scheme, graupel particles are converted to hail once wet growth occurs. Wet growth happens when graupel particles accrete an excessive amount of supercooled water, releasing a large amount of latent heat of freezing that prevents all the supercooled water on the surface from freezing. The remaining unfrozen supercooled water forms a liquid film on the particle surface. For a given particle size, whether wet growth occurs depends on the environmental temperature, humidity, and supercooled water content. It means that under the same environmental conditions, whether wet growth occurs depends on the particle size, with the critical diameter for wet growth referred to as the Schumann–Ludlam limit (SLL; Young, 1993). SLL is the diameter of a particle at which thermal equilibrium is achieved, assuming the particle surface is maintained at 0 °C. This equilibrium condition accounts for the heat balance of the sensible heat from the collected supercooled water, the sensible heat exchange between the particle surface and the surrounding air, and the latent heat of freezing and vapor diffusion at the particle surface. In the MY2 scheme, SLL is approximately calculated as follows:

$$\begin{array}{}\text{(11)}& D_{\mathrm{h}\mathrm{0}}=\mathrm{0.02}\times \left(e^{\left(\frac{-T_{\mathrm{c}}}{\mathrm{1.1}\times {\mathrm{10}}^{\mathrm{4}}\left(q_{\mathrm{c}}+q_{\mathrm{r}}\right)+\mathrm{1}}\right)}-\mathrm{1}\right),\end{array}$$

where q_c and q_r represent the mass content of supercooled cloud and rain (unit: kg m⁻³), respectively, and T_c denotes temperature (unit: °C). For a group of graupel particles, particles smaller than SLL continue to undergo dry growth, while those larger than SLL undergo wet growth and are converted to hail.

The aforementioned approach of determining whether graupel particles convert to hail based on whether wet growth occurs is physically grounded and has been adopted by many other microphysics scheme development research projects, including those by Ziegler (1985), Mansell et al. (2010), and Tsai and Chen (2020). However, this method has a drawback: although D_h0 is a function of both temperature and supercooled water content, it becomes increasingly sensitive to temperature as it approaches the 0 °C level. As shown in Khain et al. (2010b), D_h0 remains nearly constant at approximately 0.5 mm at temperatures around −2 °C, regardless of the liquid water content (Fig. 14 in their study). It means that near −2 °C, any graupel particle larger than 0.5 mm would undergo wet growth and convert to hail, no matter how much supercooled water it encounters, which is quite unreasonable.

Figure 7 illustrates the temperature and graupel size when the graupel-to-hail conversion process occurs in three heavy rainfall cases simulated with the original MY2 scheme. As shown in the figure, approximately 90 % of the graupel-to-hail conversions occur at temperatures above −0.5 °C in all three cases, and about 45 % of the graupel particles are less than 1 mm in diameter when they are converted to hail. This phenomenon of graupel particles being too easily converted to hail near the 0 °C level is attributed to the imperfect graupel-to-hail conversion parameterization mentioned above.

Figure 7 Two-dimensional frequency distribution of the location where the conversion from graupel to hail occurred as a function of temperature and the derived Schumann–Ludlam limit.

Mansell et al. (2010) identified issues with this graupel-to-hail conversion parameterization method and introduced several artificial constraints to reduce hail overforecasting. For instance, they excluded the process of graupel collecting supercooled raindrops when calculating D_h0, limited the conversion amount to no more than 1$/$10 of the graupel mass, and prevented conversions at temperatures warmer than −2 °C. Although these artificial constraints in the NSSL scheme help to reduce hail overestimation to some extent, they inevitably result in underestimation of hail occurrence and smaller hail sizes in the simulations, as indicated in Labriola et al. (Fig. 4f in 2019a, Fig. 2e in 2019b).

3.2 Results with the modified MY2 scheme

Unlike the NSSL scheme, which imposes artificial constraints lacking a physical basis, we introduce the physically grounded spongy wet growth process (described in Sect. 2.2) to address the overprediction of hail formation in the original graupel–hail conversion parameterization in the MY2 scheme. Figure 8 shows the radar reflectivity and spatial distribution of graupel and hail particles for the three heavy rainfall cases simulated with the modified MY2 scheme. Compared to the results from the original MY2 scheme shown in Sect. 3.1, the simulations provided by MY2_Mod exhibit a successful mitigation of hail overforecasting in both horizontal and vertical directions. This reduction of hail is accompanied by a decrease in radar reflectivity intensity. The simulation results provided by MY2_Mod show fewer regions with the column maximum reflectivity greater than 55 dBZ, aligning more closely with the observations. In the vertical direction, regions with reflectivity larger than 35 dBZ above the 0 °C level are primarily composed of graupel rather than hail, which is also similar to those shown in the observations.

Figure 8 Column maximum (top row) radar reflectivity and cross-sections (bottom row) of radar reflectivity simulated with the modified MY2 scheme at (a, d) 17:20 UTC on 23 June 2022, (b, e) 20:40 UTC on 17 June 2023, and (c, f) 18:00 UTC on 6 July 2023. The “×” symbols (dots) indicate the presence of hail (graupel) according to the criteria set in the paper. Gray lines represent the 0 °C level. Black lines in the top row represent the location of the cross-sections in the bottom row.

The modification of the MY2 scheme to suppress the overprediction of hail also leads to a more realistic representation of the precipitation mechanism in the simulated Meiyu systems. After the modification, the primary source of precipitation remains the accretion of raindrops, accounting for 51.5 %, 60.3 %, and 63.8 % in the three heavy rainfall cases, similar to the original MY2 scheme (Fig. 9). The contribution from hail melting is significantly reduced, dropping from 22.5 %, 17.7 %, and 17.5 % to just 0.3 %, 0.2 %, and 0.1 %, respectively. Meanwhile, the contribution from graupel melting increases substantially, rising from 20.6 %, 18.3 %, and 17.8 % to 46.2 %, 36.3 %, and 34.1 %. Although we cannot directly determine the contribution of each microphysical process to the rainwater content from the observations, the absence of hail and the dominance of graupel in regions with strong reflectivity above the 0 °C layer in the observations suggest that the modified scheme provides a more realistic representation of the precipitation processes compared to the original version. The reduction of hail is also shown in Fig. 3, where the hail percentages of all three heavy rainfall cases produced by the modified MY2 scheme are nearly zero for at least 12 h, aligning with observations. In fact, compared to the original version, the modified MY2 scheme reduces hail mass by a factor of 100 (not shown). However, it should be noted that the significant reduction in hail shown in Fig. 8 is based on the hail classification criteria outlined in the paper, and it does not imply that the hail mass produced by MY2_Mod is nearly zero. The simulations using the modified MY2 scheme do contain regions with positive hail mass mixing ratios (> 0 g kg⁻¹) that fail to meet the three specified hail identification criteria. However, these regions exhibit hail signals below the detection thresholds and are therefore negligible.

Figure 9 Same as Fig. 6 but for simulations conducted with the modified MY2 scheme.

At this stage, we have developed an improved MY2 scheme that effectively mitigates the overforecasting of hail in heavy rainfall events. However, the question remains of whether this scheme can still accurately simulate hail events. While the abovementioned regulation of wet growth in graupel particles is physically justified, it might be overly stringent, potentially hindering the simulation of hail in real hail events. Thus, we used the MY2_Mod scheme to simulate three hail events that occurred on 30 April 2021, 4 June 2022, and 25 March 2024. Figure 10 presents the observed and simulated radar reflectivity, along with the spatial distribution of graupel and hail for these events. The results show that the pattern and location of the simulated radar echoes closely match the observations. Moreover, the modified MY2 scheme successfully simulates hail in all three hailstorms, demonstrating that despite the stricter conditions for hail formation, the scheme retains its ability to produce hail.

Figure 10 Observed (first two rows) and simulated (last two rows) column maximum (first and third rows) and vertical cross-sections (second and fourth rows) of radar reflectivity for hail cases on (first column) 30 April 2021, (second column) 24 June 2022, and (third column) 25 March 2024. The “×” symbols (dots) indicate the presence of hail (graupel) according to the criteria set in the paper. Gray lines represent the 0 °C level. Black lines in the first and third rows represent the location of the cross-sections in the second and fourth rows, respectively.

3.3 Hail formation mechanism

The above analysis demonstrates that the modified MY2 scheme could distinguish between heavy rainfall and hail events, a capability absent in the original MY2 scheme and several other microphysics schemes. This modified scheme can therefore be employed to investigate why hail forms in hail events but not in heavy rainfall cases. Figure 11a illustrates the conditions where the graupel-to-hail conversion occurs in three hail cases, alongside a comparison of four associated microphysical processes. The figure shows that graupel-to-hail conversion primarily occurs when, within one integration time step, the mass of supercooled raindrops collected by graupel particles (Q_rg) exceeds the combined mass of supercooled water required for wet growth (Q_gwet) and the mass necessary to fill the internal voids of graupel particles to reach hail density (Q_fill). The conversion of graupel to hail due to sufficient collection with supercooled raindrops, rather than supercooled cloud droplets, is also justified within the MY2 scheme framework. In the MY2 scheme, the densities of graupel and hail particles are predefined: 400 kg m⁻³ for graupel and 900 kg m⁻³ for hail. To increase the density of graupel to that of hail within a single integration time step, collections with supercooled cloud droplets alone are insufficient. It can only be reached by collecting large supercooled raindrops. Additionally, the amount of supercooled water required to fill the internal voids of graupel particles (approximately 10^−6.5) is about an order of magnitude greater than that required for graupel wet growth (10^−7.5). It is this process that significantly mitigates the overly easy formation of hail in heavy rainfall events.

Figure 11 (a) Box plot of four microphysical process conversion rates (unit: log₁₀ (kg kg⁻¹ dt⁻¹)) associated with the conversion process from graupel to hail in three hail events. Dt represents the model simulation time step. Q_cg and Q_rg represent the mass of cloud droplets and raindrops collected by graupel per time step, respectively. Q_gwet stands for the critical mass of supercooled water required for graupel particles to undergo wet growth. Q_fill represents the mass of unfrozen water required to fill the internal voids of graupel particles to achieve the bulk density of pure ice. (b) Box plot of vertical velocities at locations where the conversion from graupel to hail occurs in (the first three columns) three hail cases and those at the corresponding locations in (the last three columns) three heavy rain events. (c) Comparison of Q_r-weighted vertical velocities and terminal fall velocity of raindrops at locations where the conversion from graupel to hail occurs in three hail cases. Q_r represents the mass mixing ratio of raindrops. (d) The same as (c) but for three heavy rainfall cases at corresponding locations.

In hail events, the vertical velocity at locations where graupel-to-hail conversion occurs is generally greater than 10 m s⁻¹. For example, during the hail event on 24 June 2022, the median vertical velocity reached 16 m s⁻¹ (Fig. 11b). Additionally, the regions where graupel-to-hail conversion occurs are relatively uniformly distributed in areas where radar reflectivity ranges from 45 to 55 dBZ and temperatures range from −5 to 0 °C (not shown). To investigate why graupel-to-hail conversion does not occur in heavy rainfall cases, we randomly sampled regions in the simulation results of three heavy rainfall events with radar reflectivity between 45 and 55 dBZ, temperatures between −5 and 0 °C, vertical velocity exceeding 0 m s⁻¹ , and the presence of supercooled raindrops. It was found that, in heavy rainfall cases, the vertical velocity in regions corresponding to those where graupel-to-hail conversion occurs in hail events is significantly lower, generally less than 10 m s⁻¹. Additionally, in heavy rainfall events, the vertical velocity does not significantly exceed the terminal fall velocity of raindrops, limiting the upward transport of raindrops to graupel regions and mitigating the hail formation process. In contrast, in hail events the vertical velocity in graupel-to-hail conversion regions generally exceeds the terminal fall velocity of raindrops, facilitating the upward transport of raindrops to graupel particles, thereby initiating the graupel-to-hail conversion process.

Why is the vertical velocity near the 0 °C isotherm in heavy rainfall events insufficient to transport raindrops to higher levels? Given that the dynamic vertical pressure gradient force (about 0.03 m s⁻²) is about an order of magnitude smaller than the buoyancy force (about 0.2 m s⁻²) in these six cases (not shown), our study primarily focuses on the dynamical effects of buoyancy. Buoyancy is calculated as

$$\begin{array}{}\text{(12)}& B=g\left[{\displaystyle \frac{\mathit{\theta }-\overline{\mathit{\theta }}}{\mathit{\theta }}}+\mathrm{0.61}\left(q_{\mathit{\nu }}-\overline{q_{\mathit{\nu }}}\right)-q_{\mathrm{hydro}}\right],\end{array}$$

where θ is the potential temperature and q_hydro the mixing ratio of multiple hydrometeors. The overbar in Eq. (12) represents the base state defined as the areal mean in the 100 km × 100 km region enclosed by the gray dashed lines in insets in the top right corner of each panel shown in Fig. 12. The explanation for why the gravitational force of hydrometeors (the last term in Eq. 12) also affects buoyancy is given below: given the short time required for hydrometeors to reach their terminal fall velocity during descent, microphysics schemes assume that hydrometeors descend at a constant terminal velocity, where air resistance is balanced by gravitational force. According to Newton's second law, the drag force exerted by the hydrometeors on the air (i.e., negative buoyancy) equals the air resistance they encounter, which is also equal to their gravitational force.

Figure 12 presents vertical profiles of vertical velocity and buoyancy at representative moments and regions for both heavy rainfall and hail events. The figure demonstrates that the main updraft regions are closely aligned with regions of high buoyancy resulting from potential temperature differences, indicating that upward motion in both hail and heavy rainfall events is primarily driven by buoyancy induced by these temperature differences. In the hailstorm cases, particularly in Hailstorm1, a negative buoyancy zone formed by the cold pool is present below 2 km altitude. The quasi-upright updraft structure indicates a dynamic equilibrium between the cold pool and environmental shear, consistent with the Rotunno–Klemp–Weisman (RKW) theory (Weisman and Rotunno, 2004), which favored the sustained severe convection. In hail events, buoyancy within the main updraft regions typically reaches 0.16 m s⁻², exceeding the average value of about 0.10 m s⁻² in heavy rainfall events. Consequently, the maximum updraft velocity in hail events (35–40 m s⁻¹) is significantly greater than that in heavy rainfall events (15–20 m s⁻¹). However, as mentioned before, the vertical velocity near the 0 °C isotherm is more crucial for graupel-to-hail conversion, and thus it is necessary to focus on buoyancy and vertical velocity near the 0 °C level.

Figure 12 Vertical cross-sections of buoyancy (shaded, unit: m s⁻²) caused by potential temperature differences, vertical velocity (black solid lines), and negative buoyancy due to cloud droplets (violet), raindrops (blue), snow (grape), and graupel (orange) in three (a, b, c) heavy rainfall and (d, e, f) hail events. The values of negative buoyancy caused by different hydrometeors are contoured at 0.01, 0.02, 0.04, 0.06, and 0.08 m s⁻². Vertical velocity contours range from 5 to 40 m s⁻¹ at intervals of 5 m s⁻¹. The gray dashed line indicates the 0 °C isotherm. Insets in the top right corner of each panel show the column maximum radar reflectivity at the corresponding time and the location of the vertical cross-section (black line). The region enclosed by the gray dashed lines represents the area used to calculate the environmental average potential temperature.

In the three heavy rainfall events analyzed in this study, the 0 °C isotherm is generally located between 5–5.5 km in altitude. The high mass mixing ratio of raindrops between the 0 °C level and the surface generates negative buoyancy of 0.04–0.08 m s⁻², which strongly inhibits the updraft below the 0 °C isotherm in heavy rainfall events. In contrast, the situation differs in hail events. For example, in the squall line event on 30 April 2021, precipitation occurs not in the main updraft region but rather in the stratiform cloud behind it, preventing the rainwater from suppressing the updraft. In the hail cases on 24 June and 25 March 2024, while precipitation does occur in the main updraft region, and the negative buoyancy caused by rain drag reaches 0.04 m s⁻², it does not completely cover the region beneath the main updraft as it does in heavy rainfall events. In hail events, precipitation in the main updraft region typically occurs at the updraft's edges (Fig. 12e), or the updraft itself is slightly tilted (Fig. 12f), preventing the drag force of rainwater from fully weakening the updraft near the 0 °C isotherm.

4 Conclusions

Many microphysics schemes that predefine separate ice-phase hydrometeor categories, such as the MY2, NSSL, NTU, and SBM schemes, consider the occurrence of wet growth of graupel particles to be a prerequisite for their conversion to hail. While this approach is physically grounded, using the critical amount of supercooled water required to raise the graupel surface temperature to 0 °C as the sole trigger for the graupel-to-hail conversion process is insufficient, since less supercooled water is required for the occurrence of wet growth as graupel particles approach the 0 °C level. Consequently, this parameterization method implies that graupel particles are more likely to be converted to hailstones near the 0 °C level. In this study, we used the WRF model configured with the original MY2 scheme to simulate three heavy rainfall cases within the Meiyu front systems in eastern China, where hail was not detected by a hydrometeor classification algorithm based on dual-polarization radar data. It is found that the original MY2 scheme indeed produced a large amount of spurious hail for the reasons mentioned above. Although these overpredicted hailstones do not reach the surface, the issue of hail overprediction cannot be considered trivial. The overpredicted hailstones produce reflectivity values higher than those observed and contribute excessively to rainfall. This, in turn, complicates our understanding of the radar echo structure, spatial distribution of hydrometeors, and precipitation formation mechanisms in the actual precipitation system. We argue that the hail overprediction shown in the paper is not confined to the MY2 scheme. Any scheme utilizing the aforementioned graupel-to-hail conversion parameterization method without proper constraints is likely to encounter similar issues of hail overprediction.

By incorporating the spongy wet growth process of graupel particles, we stipulate that graupel-to-hail conversion occurs only when the amount of supercooled cloud water and rainwater collected by graupel particles within one integration time step exceeds the amount required to maintain their surface temperature at 0 °C and to fill the internal voids of the graupel particle so that its density reaches the density of pure ice. The modified MY2 scheme greatly mitigated the overprediction of hail. Furthermore, through simulations of three typical hail events in eastern China, the modified MY2 scheme still successfully simulated the presence of hail in hail events, demonstrating its ability to distinguish between heavy rainfall and hail processes – an ability that the original MY2 scheme and many other schemes lacked.

We further compared the heavy rainfall events and hail events simulated by the modified MY2 scheme to analyze why the former does not lead to hail formation, while the latter does. It is found that hail formation primarily occurs when the amount of supercooled rainwater collected by graupel particles exceeds the amount required to fill the internal voids of the graupel particles. Therefore, the transport of rainwater across the 0 °C level into the subfreezing region is crucial for hail formation. Compared to the three heavy rainfall cases, the three hail cases exhibit stronger updrafts driven by greater positive buoyancy resulting from the temperature difference between the convective updraft and the surrounding environment. Additionally, in the hail cases, the primary updraft region is spatially separated from precipitation, and thus the rainwater drag force could not completely inhibit the updraft. In contrast, in the heavy rainfall cases, precipitation directly occurs beneath the convective updraft. Therefore, the negative buoyancy from rainwater drag weakens the updraft under the 0 °C layer, causing the updraft velocity to generally be less than the terminal fall velocity of raindrops, preventing them from being transported to the subfreezing layer to be collected by graupel particles and forming hail.

5 Discussion

Our proposed method of introducing spongy wet growth to improve the graupel-to-hail conversion process still requires further investigation in several respects. For instance, the model integration time step may exert an influence on the conversion occurrence. Larger dt values increase the likelihood that the supercooled water content accreted by graupel particles exceeds the amount required to fill their internal voids, thereby artificially enhancing hail formation. Although our current experiments employ a sufficiently small dt (2 s), comprehensive sensitivity tests are needed in the future to provide optimal dt thresholds. Furthermore, the predefined graupel and hail density in the MY2 scheme physically constrains hail conversion by requiring graupel particles to accrete sufficient supercooled water within a single time step. This inherent configuration may systematically suppress hail production. Future implementation of this approach within microphysical schemes which could predict the density of graupel and hail particles, such as the prognostic graupel density MY2 (Milbrandt and Morrison, 2013) and the NSSL scheme, could potentially mitigate this limitation. Additionally, the absence of high-density ice particles above the freezing level, as shown in the P3 scheme (Morrison and Milbrandt, 2015; Milbrandt and Morrison, 2016; Milbrandt et al., 2021) by Labriola et al. (2019b, Fig. 14a in their study) and the Thompson_Hail scheme (Jensen et al., 2023) by Li et al. (2024, Fig. 8f in their study), could potentially be improved by introducing the spongy wet growth process to increase the density of riming particles.

Previous studies have primarily relied on real hail case studies to evaluate HMSs, focusing on whether these schemes can accurately reproduce radar structures, simulate hail reaching the ground, match observed surface hail locations, and estimate surface hail size. These studies are crucial as they highlight the sources of simulation biases in hail prediction and have led to improvements in microphysics schemes. In this paper, we propose an alternative evaluation approach: assessing hail microphysics schemes in heavy rainfall events where no hail is present – not only is there no hail reaching the ground, but there is also no detectable hail aloft, as determined by dual-polarization radar. The goal is to test whether the scheme correctly refrains from simulating hail in such cases. We argue that a robust hail microphysics scheme should accurately simulate hail in real hail events, while refraining from simulating it in storms where no detectable hail exists. This approach offers valuable insights for the development of microphysics schemes and for advancing our understanding of hail formation. We recommend using these hail-free cases to evaluate other HMSs, including those that predefine ice-phase particle categories (e.g., NSSL) and allow particle properties to evolve freely (e.g., different configurations of the P3 scheme). This would help assess their ability to accurately simulate both hail and non-hail processes and identify potential areas for further model improvements.