The squall-line event on 20 May 2011, during the Midlatitude Continental Convective Clouds (MC3E) field campaign has been simulated by three bin (spectral) microphysics schemes coupled into the Weather Research and Forecasting (WRF) Model. Semi-idealized three-dimensional simulations driven by temperature and moisture profiles acquired by a radiosonde released in the preconvection environment at 1200 UTC in Morris, Oklahoma, show that each scheme produced a squall line with features broadly consistent with the observed storm characteristics. However, substantial differences in the details of the simulated dynamic and thermodynamic structure are evident. These differences are attributed to different algorithms and numerical representations of microphysical processes, assumptions of the hydrometeor processes and properties, especially ice particle mass, density, and terminal velocity relationships with size, and the resulting interactions between the microphysics, cold pool, and dynamics. This study shows that different bin microphysics schemes, designed to be conceptually more realistic and thus arguably more accurate than bulk microphysics schemes, still simulate a wide spread of microphysical, thermodynamic, and dynamic characteristics of a squall line, qualitatively similar to the spread of squall-line characteristics using various bulk schemes. Future work may focus on improving the representation of ice particle properties in bin schemes to reduce this uncertainty and using the similar assumptions for all schemes to isolate the impact of physics from numerics.
A spectral or bin microphysics scheme predicts hydrometeor concentrations in discrete mass or diameter bins, while a bulk scheme predicts one or more moments of bulk quantities of the size or mass distribution. It is often believed that spectral or bin microphysics schemes, by having more degrees of freedom in representing hydrometeor size distributions and calculating microphysical processes, provide more realistic simulations and reduced uncertainty compared with bulk microphysics (Khain et al. 2015). This has been the primary justification for their use, given the much greater computational cost compared to bulk schemes. As such, bin schemes have been widely used to develop new bulk parameterizations or improve existing ones by providing parameterized relationships of microphysical processes based on bin simulation results (e.g., Khairoutdinov and Kogan 2000; Thompson et al. 2004; Seifert and Beheng 2006; Seifert et al. 2006; Morrison and Grabowski 2007; Saleeby and Cotton 2008; Li et al. 2009a,b; Milbrandt and McTaggart-Cowan 2010; Shipway and Hill 2012; Kogan and Belochitski 2012; Wang et al. 2013; Khain et al. 2015; Igel and van den Heever 2017a,b).
The earliest bin schemes calculated drop size distribution functions defined on a mass grid (e.g., Berry 1967; Berry and Reinhardt 1974; Bleck 1970; Soong 1974). These functions obeyed the normalization condition, namely, the integral over the particle spectrum is equal to the total particle number concentration. Schemes were improved by the development of more sophisticated numerical techniques, including the prediction of two moments in each bin (Tzivion et al. 1987; Feingold et al. 1988; Tzivion et al. 1999) and solving particle size distribution (PSD) equations over mass grids conserving multiple moments of the distribution (Khain et al. 2004). Over the past few decades, bin schemes have been extended to include both liquid- and ice-phase microphysics (e.g., Takahashi 1976; Hall 1980; Khvorostyanov et al. 1989; Khain et al. 1993; Khain and Sednev 1995, 1996; Geresdi 1998; Geresdi et al. 2014; Reisin et al. 1996; Khain et al. 2004; Lebo and Seinfeld 2011; Lebo et al. 2012). Ice processes are considerably more challenging to represent than liquid due to the wide variety of ice crystal habits, precipitating ice particle shapes and characteristics (Hashino and Tripoli 2007). While bin schemes that include ice hydrometeors can model the evolution of ice particle spectra with more confidence, they still have uncertainties in assuming ice particles properties such as density and fall speed.
vanZanten et al. (2011) recently compared different bin schemes within fully dynamical cloud models. In the simple context of liquid-only microphysics, there was a large spread for quantities such as surface precipitation rate using different bin schemes, even greater than the differences between simulations using bulk microphysics. This raises important questions about the robustness of using bin schemes to improve bulk schemes. While to our knowledge few detailed intercomparison studies between different mixed-phase bin microphysics schemes have been conducted, one might anticipate an even greater spread given the large uncertainties pertaining to ice microphysics; however, Ovchinnikov et al. (2014) found reasonable agreement between results of the System for Atmospheric Modeling spectral (bin) microphysics (SAM-SBM) (with the microphysics based on Khain et al. 2004) and Distributed Hydrodynamic Aerosol and Radiative Modeling Application (DHARMA-bin) simulations for an Arctic stratiform cloud case. It is important to point out that both vanZanten et al. (2011) and Ovchinnikov et al. (2014) compared schemes in different dynamical models, and it is difficult to ascribe their spread solely to the representation of cloud microphysics.
This study compares three different bin microphysics schemes in simulations of moist deep convection using the same dynamical model. We simulate a squall-line mesoscale convective system, which occurred on 20 May 2011, in the central United States during the Midlatitude Continental Convective Clouds Experiment (MC3E; Jensen et al. 2016; Wu and McFarquhar 2016). A squall-line case was chosen for several reasons. First, mature squall lines contain a wide variety of ice hydrometeor types (unrimed small ice, large aggregates, graupel, etc.) evolving through various growth mechanisms, and hence serve as an ideal test bed for microphysical schemes. Second, previous studies have shown large sensitivity of squall lines to the representation of microphysics because of tight coupling between the microphysics and storm dynamics (e.g., Fovell and Ogura 1989; McCumber et al. 1991; Ferrier et al. 1995; Liu et al. 1997; Morrison et al. 2009, 2012). Finally, there have been several previous intercomparison studies of bulk schemes for squall lines (Fovell and Ogura 1989; McCumber et al. 1991; Morrison et al. 2015; Fan et al. 2015, 2017), providing a context for the current bin scheme intercomparison.
Previous studies have demonstrated a large spread of squall-line simulations using different bulk schemes (Fovell and Ogura 1989; McCumber et al. 1991; Morrison et al. 2015; Fan et al. 2015). This spread has been attributed to the impact of evaporation and melting rates affecting the evolution of cold pools (Morrison et al. 2012, 2015). Rotunno et al. (1988, hereafter RKW88) provided a context for the impact of changes in cold pool properties on the storm dynamics, the so-called Rotunno–Klemp–Weisman (RKW) theory (RKW88). This theory has been revised and updated in several recent studies (Weisman and Rotunno 2004; Bryan et al. 2006; Bryan and Rotunno 2014). The RKW theory posits that an “optimal state” occurs when the vorticity associated with the cold pool is balanced by that associated with the low-level environmental shear, leading to the generation of upright, strong convective updrafts rather than tilted, relatively weak updrafts (Parker 2010). Thus, changes in the microphysics affecting latent cooling drive differences in cold pools and hence dynamics through the mechanisms outlined by RKW theory.
In this study, the macrophysical, thermodynamic, and dynamical characteristics of squall-line simulations using three bin schemes are compared. Results are discussed in the context of detailed in situ and remote sensing observations, and analyzed within the framework of RKW theory. The primary goal is to characterize the spread of simulated storm characteristics using various bin schemes in their default configurations and to understand key drivers of these differences. Implications and recommendations for the use of bin schemes to test and improve bulk schemes are also discussed.
The remainder of the paper is organized as follows. Section 2 briefly describes the three bin schemes. Section 3 presents the experimental design. Observations and model results are compared in section 4. The cold pool dynamics and the microphysics–dynamics interactions are analyzed in detail within the framework of RKW theory in section 5. Conclusions and recommendations for future work are provided in section 6.
2. Model description
The three bin microphysics schemes applied in this study have been implemented into the Weather Research and Forecasting (WRF) Model, version 3.4.1 (Skamarock et al. 2008). The key features of each scheme related to the current work are briefly described in this section. Some details about each scheme such as the mass–size relation, terminal velocity, and density function are documented in Table 1. Interested readers are recommended to look at the details of each scheme in the related references.
a. CNNB scheme
The Caltech–NCAR–NOAA Bin (CNNB) scheme (Lebo and Seinfeld 2011; Lebo et al. 2012) is a detailed mixed-phase two-moment bin microphysics scheme that applies the method of moments (Tzivion et al. 1987) and combines the formulations of Tzivion et al. (1987, 1999), Feingold et al. (1988), Reisin et al. (1996), and Khain et al. (2004) (see details in Lebo and Seinfeld 2011). The CNNB scheme applies a mass-doubling bin approach to simulate evolution of the number and the mass mixing ratios of water, ice, snow, and graupel/hail in 36 bins. Aerosols are treated prognostically in the model, being advected and mixed by the flow as well as being removed and regenerated upon droplet activation and complete evaporation, respectively. The particles are distributed following a log-normal size distribution with a geometric mean diameter of 0.1 μm and a log standard deviation of 1.8. Activation is predicted using Kohler theory to determine the critical supersaturation of aerosol particles in each bin. A time-dependent melting model based on Rasmussen and Heymsfield (1987) has recently been implemented for ice, snow, and graupel. Since shedding of water from melting ice particles does not occur until ice particle size is larger than several millimeters, in most cases the melted water stays with the melting ice particle during sedimentation. New prognostic variables for the melted water mixing ratio were added for each ice, snow, and graupel/hail bin to track the partial melting. For ice nucleation, homogeneous and immersion freezing are parameterized following Bigg (1953) and Vali (1975); condensation and deposition follow Meyers et al. (1992) and Harrington et al. (1995), while the secondary ice formation mechanism from rime splintering (Mossop and Hallett 1974) is implemented following prior works (e.g., Beheng 1987). Hydrometeors are assumed spherical and the densities of snow and graupel/hail particles are the highest among all of the schemes tested. CNNB contains more than 400 scalar variables and uses a small microphysics time step for all microphysical processes. The computational cost is about 100–400 times more than a typical bulk scheme simulation.
b. FSBM scheme
The Fast Spectral-Bin Microphysics (FSBM) is a simplified version of the Hebrew University Cloud Model (HUCM) described by Khain and Sednev (1996), Khain et al. (2004), and Fan et al. (2012). Both the full HUCM SBM and FSBM are included within the WRF official release and have been widely used by the community. The choice of FSBM for this study is due to its much lower computing cost compared to full HUCM SBM and its lower technical barrier to implement necessary changes specific for this intercomparison. The original microphysics parameterization (i.e., the full version of HUCM SBM) solves a system of kinetic equations for the size distribution functions of water drops, ice crystals (plate, columnar, and branch types), snow/aggregates, graupel and hail/frozen drops, as well as cloud condensation nuclei (CCN). The PSD is calculated on a finite-difference mass grid containing 33 or 43 mass-doubling bins. Unlike the other two schemes utilizing the method of moments (Tzivion et al. 1987), FSBM applies a bin-resolved moment-conserving technique (Khain et al. 2004, 2015) to solve PSD equations for each size bin and both number and mass are predicted at the same time. The FSBM contains four size distributions (CCN, water drops, ice/snow, and graupel), derived from the full HUCM by ignoring the three ice crystal and hail categories. It calculates ice and snow using a single size distribution divided at a diameter of 300 μm. The CCN size distribution is prognostic with sinks and sources, which include CCN sources at boundaries, advection, and nucleation for the version in the WRF Model. Scavenging of CCN by precipitation and CCN regeneration upon evaporation of liquid drops are not considered. Ice nucleation parameterizations for FSBM are the same as in the original HUCM documented in Khain et al. (2004) (i.e., Bigg 1953) for both heterogeneous and homogenous freezing of droplets, Meyers et al. (1992) for condensation and deposition nucleation with a maximum ice supersaturation limit at 50%, and Mossop and Hallett (1974) for the secondary ice formation. There are several versions of the full SBM with different bin numbers and levels of complexity, as well as FSBM, which incorporate different assumptions of hydrometeor properties for high-density ice particles for studying different types of clouds (Fan et al. 2015, 2017). The code used in this study is the default version in the WRFV3.6.1 release from 2014 (which we back ported to WRFV3.4.1). The default version is chosen for comparing since it is the version that can be publicly accessed and is used extensively. FSBM considers the shape of the hydrometeors implicitly when calculating fall velocities and collision kernels. For example, the aspect ratio of 0.8 is assumed for snow particles whereas a ratio of 1 (spherical shape) is assumed for all other hydrometeor species for collisional cross-section calculations. The densities of snow and graupel particles are the lowest among all the schemes tested. FSBM has the fewest scalar variables among three schemes (~150) and costs about 20–80 times more than a bulk scheme simulation but about 4 times less than the CCNB and UPNB schemes.
c. UPNB scheme
The University of Pécs and NCAR Bin (UPNB) scheme (Geresdi 1998; Rasmussen et al. 2002; Xue et al. 2010, 2012; Geresdi et al. 2014) is a detailed mixed-phase two-moment bin microphysics scheme applying the method of moments (Tzivion et al. 1987, 1999; Feingold et al. 1988; Reisin et al. 1996) and mass-doubling approach to simulate the evolution of size distributions of water, pristine ice, snow, and graupel with 36 bins in each category. The UPNB scheme tracks the rime fraction of snow (the ratio of the rimed water mass to the total snow particle mass) to smoothly convert heavily rimed snow to graupel (Sarkadi et al. 2016), which is a unique feature among the three tested schemes. A detailed melting algorithm based on Rasmussen and Heymsfield (1987) was implemented for snow and graupel (Geresdi et al. 2014; Sarkadi et al. 2016). This scheme has the following prognostic variables to represent the aforementioned processes: the mass of rimed water on snowflakes, the number concentration of the snow particles that have collected at least one drop, and the mass of melted water on melting snow and on melting graupel particles in each bin. A size-resolved aerosol activation and regeneration scheme based on the Köhler theory and the hygroscopicity concept has been developed for the UPNB scheme (Xue et al. 2010, 2012; Muhlbauer et al. 2010; Chen et al. 2011; Teller et al. 2012). However, in the current study, a cumulative CCN spectrum as a function of supersaturation respect to water was used to simulate the cloud droplet activation, which is similar to the CCN activation method in FSBM. The homogeneous freezing process is described by Bigg (1953), the condensation freezing and deposition nucleation processes are parameterized by the modified Cooper curve (Cooper 1986), and secondary ice formation follows Mossop and Hallett (1974). All ice-phase particles assume a nonspherical shape in UPNB. Similar to CNNB, UPNB also contains more than 400 scalar variables. The computational cost is similar to CNNB and about 4 times more than FSBM.
In summary, these three schemes all include key microphysical processes leading to the formation of liquid and ice cloud and precipitation particles. However, there are substantial differences in the parameterized physics and in the numerical methods for solving the governing equations. These include the following: different representations of the PSD, method of moments versus a moment-conserving flux-based method for representing PSD evolution, different treatments of ice nucleation and melting, different mass ranges, and different mass-size (density) and terminal velocity–size relationships. The differences in mass–diameter relationships, terminal velocities, and ice nuclei (IN) concentrations for the three schemes are shown in Fig. 1. CNNB assumes the smallest sizes and the highest fall speeds for snow and graupel due to their assumed spherical shapes and high densities. FSBM has the lowest fall speeds for snow and graupel particles while UPNB has the largest maximum particle sizes; UPNB’s fall speeds for snow and graupel are between CNNB and FSBM. Because there are fewer bins and larger mass in the first bin, the maximum mass of hydrometeors in FSBM is about one-fifth of the other two schemes. UPNB also has weaker temperature dependence of the heterogeneous ice nucleation than the other two schemes (cf. Fig. 1c).
The differences in model results simulated by these schemes reflect the combined effects of the hydrometeor properties, parameterized physics, and numerics. To test the impact of each individual aspect, one can run the case with just one scheme using different assumptions or run the case with multiple schemes using the same assumptions. For example, in order to test the sensitivity of the simulated squall line to the numerics, one can assume the FSBM mass–size and fall speed relations for CNNB and UPNB. Such experiments are beyond the scope of this study. Future work will focus on separating the impact of each individual aspect of the schemes.
3. Experimental design
A three-dimensional (3D) semi-idealized model setup1 similar to Bryan and Morrison (2012) and Morrison et al. (2012) is used in this study. Open boundary conditions in the across-line (x) direction and periodic boundary conditions in the along-line (y) direction are specified in a 3D domain with a size of 612 km (x) by 122 km (y) by 25 km (z). The horizontal and vertical grid spacings are 1 km and approximately 250 m, respectively. The dynamic time step is set to 3 s. Rayleigh damping is applied to the upper 5 km.
The main idea is to investigate sensitivity of the simulated squall line to different bin microphysics schemes in the context of observations but in a simplified framework to facilitate the scheme comparison and focus on microphysical–dynamical interactions, similar to previous studies using bulk microphysics (Bryan and Morrison 2012; Morrison et al. 2012, 2015). Thus, we use an idealized setup and turn off the planetary boundary layer (PBL), land surface, and radiation schemes. All simulations are integrated for 6 h with the first 2 h as the spinup time. We are aware of the sensitivity of simulated squall-line properties to the grid spacing (cf. Bryan and Morrison 2012), but the computational expense of bin microphysics prevented us from running additional simulations at finer resolutions.
The simulations were driven by an environmental sounding based on observations from Morris, Oklahoma, at 1200 UTC 20 May 2011 (Fig. 2). Therefore, the diurnal temperature variation during the 6 h of squall-line evolution is not considered in this work. Figure 2b shows the smoothed thermodynamic profiles and the idealized wind profile of the model sounding. This preconvection sounding features relatively stable low-level air below 900 m, very moist air below 3 km, and relatively steep lapse rates between 4 and 6 km with CAPE of 2200 J kg−1 from 900 m up. The winds are mostly from the south-southwest below 4 km and the shear perpendicular to the squall line is about 0.003 s−1 below 4 km (not shown). This is approximated by using a unidirectional shear (in the x direction) of 0.003 s−1 between 0 and 4 km (U from −12 to 0 m s−1) (Fig. 2b). For simplicity, there is no wind above 4 km in the simulation, which is not the case in the real storm.
Unlike Morrison et al. (2012), in which a cold pool was prescribed to initiate the convection, a time-varying low-level horizontal momentum convergence is applied for the convection initiation in this study (Loftus et al. 2008; Morrison et al. 2015). The horizontal convergence is fixed at the center of the x direction with a maximum acceleration of 0.1 m s−2 uniformly distributed along the y direction at the lowest level. The acceleration is radially reduced to 0 from the center to 10 km away in ±x and +z directions. This across-line convergence pattern is applied to the simulations during the first 55 min and then linearly decreased to zero at 1 h. Since the cold pool strength plays a key role in determining the squall-line structure, we view the convergence method for convective initiation as an improvement over the cold pool insertion method since it avoids imposing artificial cold pool properties in the simulations.
For CNNB, a background aerosol concentration of 320 cm−3 with a single mode lognormal distribution is used to initiate cloud droplet activation based on the CCN concentration at 0.4% supersaturation and cloud droplet observations of this case (Fan et al. 2015). The aerosol composition is assumed to be ammonium sulfate. For FSBM and UPNB, CCN spectra as functions of supersaturation with respect to water are prescribed to generate similar cloud droplet concentrations to CNNB of about 300 cm−3. We note that this value is low for a typical continental storm, especially in the convective region, which may produce uncertainties in simulating the convective properties of this case.
4. Model–observation comparisons
In this section, the squall-line storm structure, vertical velocity, precipitation, and cold pool properties simulated by the three bin schemes are compared to the observed NEXRAD S-band radar reflectivity [3D clutter-removed Ze data (Fan et al. 2017)], the vertical velocity derived from both the CSAPR C-band and the X-band radars (Fan et al. 2015; Collis et al. 2013) and the mesonet observations (Oklahoma Mesonet 2015). As stated previously, the main purpose of this study is not to reproduce the observed case but to investigate the sensitivity of simulated squall-line structure to different bin schemes. Thus, the model–observation comparison does not serve to validate the schemes but provides a context for the simulation differences. While previous studies showed that semi-idealized simulations and observations can be reasonably compared (Bryan and Morrison 2012; Morrison et al. 2012, 2015), only the main features of the storm are captured by these simulations since many physical processes (large-scale flow patterns, PBL, land surface and radiation) are omitted.
a. S-band radar reflectivity comparisons
The 3D gridded NEXRAD reflectivity data at 1200 UTC within the red box in Fig. 2a have been processed to generate the line-averaged dBZ vertical cross section (average of dBZ field in the along-line direction) and the horizontal cross section at 2 km in Figs. 3a1 and 3b1. This region was subjectively selected to cover the well-developed linear convective line and widespread stratiform region in a similar along-line dimension to the model domain so that the along-line-averaged fields are comparable to the model results. Because the observed squall line at 1200 UTC was well into its mature stage, the modeled squall lines at the end of the simulation (hour 6), when they are also well into their mature phase, are plotted in Figs. 3a2–4 and 3b2–4 for CNNB, FSBM and UPNB, respectively. For the model reflectivity, an algorithm similar to Sarkadi et al. (2016) was applied. This algorithm, based on Blahak (2007) and Smith (1984), explicitly considers the dependence of reflectivity on melted fraction and on the size of the melting particles. Since FSBM does not predict the melt fraction, this effect is neglected in the FSBM reflectivity calculations. Hydrometeor size distributions simulated by the bin schemes, the specific radar wavelength (10 cm in this case), temperature, and phase-dependent refractivity are also considered (Ray 1972; Smith 1984; Mätzler 1998). Details about this radar reflectivity algorithm for bin schemes can be found in Sarkadi et al. (2016).
A subsection of the model domain (x from100 to 550 km and z from 0 to 15 km) is shown in Figs. 3a2–4 to compare with the observed reflectivity (Fig. 3a1) in the same dimension that covers the entire extent of the simulated squall lines. Each scheme simulates a squall line with a strong convective leading edge, a transition zone and a trailing stratiform region, which is qualitatively similar to the observed reflectivity. In the convective core region, high reflectivity (>30 dBZ) was mainly observed below 4 km (about the 0°C isotherm based on sounding in Fig. 2b), indicating enhanced reflectivity by melting particles in NEXRAD data, while all of the schemes simulate high reflectivity above this altitude. On the other hand, the observed high reflectivity below the melting layer in the stratiform region is not captured by the models. Although CNNB simulates a deep melting layer extending from 4 km to the ground, the lack of hydrometeors and their relatively smaller sizes in the stratiform region above 4 km results in weak reflectivity in the corresponding melting layer. FSBM simulates a broader stratiform region with higher reflectivity than the other two schemes, closer to observations. Since FSBM does not explicitly simulate the melted water coating on ice particles, the reflectivity is not enhanced within the melting layer. A well-defined bright band between 3 and 4 km with reflectivity comparable to the observations is simulated by UPNB. However, most ice particles are completely melted at 3 km with weaker reflectivities below compared to the observations due to strong evaporation.
The 120-km-wide horizontal cross section at 2 km AGL of the southwestern portion of the observed data are plotted in Fig. 3b1 to compare with the corresponding model results. At this level, reflectivities as high as 55 dBZ were observed in the convective leading edge and even higher reflectivity is simulated by CNNB due to large, water-coated, partially melted ice particles. Large spatial variability is seen in both the observed and simulated reflectivity fields. Averaging along the line (i.e., in the y direction) smooths the field and reduces the maximum reflectivity in the vertical cross-section plots. This is especially true for the observations because the observed convective line was not as linear as its simulated counterparts. FSBM simulates the widest stratiform region with weakest reflectivity in the convective region among all schemes. The features are mainly attributed to the low-density and slowly falling graupel assumed in FSBM that will be discussed further. The maximum reflectivity in the stratiform region is around 30–35 dBZ, which is similar to observations. However, the area of reflectivity greater than ~30 dBZ and the overall width of the stratiform region are underestimated by all schemes. It is possible that the simulated stratiform region could grow wider given a longer simulation time. However, longer simulations are not possible with this domain configuration because the storm moves too close to the right lateral boundary.
In summary, the S-band reflectivity comparison shows that the bin schemes are able to simulate a squall-line system qualitatively similar to observations. However, there are significant differences in the general structure of the simulated system. Explicit treatment of partially melted ice particles in the scheme seems necessary to capture the observed high reflectivity in the melting zone of the stratiform region.
b. Vertical velocity comparisons
On 20 May 2011, both the C-band and the X-band radars worked in the SGP inner grid (as indicated by the blue box in Fig. 2a) before 1100 UTC. From 0800 to 1100 UTC, the northern part of the squall line passed through the area from southwest to northeast. The vertical velocity has been retrieved by combining both the C-band and the X-band radar signals using the simplified 3D variational technique following Collis et al. (2013). Only the updraft cores (vertical velocity greater than about 2 m s−1) were retrieved with confidence. Contoured frequency by altitude diagrams (CFADs) of vertical velocity greater than 2 m s−1 are plotted in Fig. 4a in logarithmic scale for the radar derived data between 0800 and 1100 UTC within the inner grid. The top at 12 km was chosen to exclude low-quality data above this level. The gridded data of radar-retrieved vertical velocity has a resolution of 500 m in both the horizontal and vertical. For the model results, vertical velocities greater than 2 m s−1 throughout the entire simulation and over the whole domain were used to construct the CFADs in Figs. 4b–d to cover the various conditions observed by the radar. Comparisons of processed radar and model data at the same resolution (1 km in horizontal and 500 m in vertical) show similar results to the raw data. Therefore, only results using the original data resolutions are presented.
Overall, the CFADs from the model simulations are similar to each other, although they differ from the retrievals partially because the bowing of the squall line and many shallow convective clouds in the real system were not simulated. None of the bin simulations capture the wide spread of the updraft distribution below 4 km. At 4-km height, FSBM and UPNB simulate slightly wider distributions than CNNB. Radar data show that vertical velocities greater than 12 m s−1 persist above 1 km, while model results only show such strong updrafts at higher altitudes. The inability of the model to capture large vertical velocities at low levels is likely because small-scale shallow convection cannot be resolved with 1-km horizontal grid spacing. Another possible reason is that the land surface, boundary layer, and radiation processes, which might affect the low-level updrafts are neglected in our semi-idealized setup. All schemes simulate peak updraft vertical velocities between 7 and 10 km, which are close to the radar estimates. The most frequent updraft velocity simulated within this layer is around 12 m s−1. This is larger than the radar values, which indicates that the models may overestimate updraft velocity in the convective cores consistent with results of Varble et al. (2014). Such an overestimate may result from unresolved turbulent mixing and lack of entrainment in updrafts at 1-km grid spacing (Bryan and Morrison 2012; Kurowski et al. 2014; Lebo and Morrison 2015), or other factors. Although the upper-tropospheric CFADs for the observations and simulations feature similar widths, the shape of the distribution, similar in all models, is substantially different than the observations.
Each scheme simulates a mode below 2 km with updrafts between 2 and 4 m s−1 for CNNB and between 2 and 5 m s−1 for FSBM and UPNB. A similar pattern is also evident from the radar data. This feature is likely the result of the mechanical lifting of air in front of the leading edge of the cold pool. The narrower distribution of this mechanical updraft mode using CNNB relative to the other two simulations suggests that the cold pool strength is weaker in this simulation (documented later in the paper). A second mode in the range of 5–8 m s−1 extending from 2 to about 3 km is simulated by each scheme. This updraft extension is associated with the latent heat release of the initial condensation resulting from the cold pool lifting. Another common feature simulated by the schemes is a mode at 6 km with vertical velocity around 12 m s−1. The velocity frequency maximum increases from 3 km and upward as a result of the convective energy release associated with the steep lapse rates between 4 and 6 km (Fig. 2b). The velocity at the frequency maximum decreases above this altitude as a consequence of more stable air aloft (the tropopause height is around 12 km). Overall, the CFAD differences between the observations and model simulations may come from various reasons such as the unrepresentativeness of the initial sounding and the relatively coarse grid spacing, although the fact that the CFADs are similar in all model simulations is encouraging. This similarity implies that differences between the schemes do not have a large impact on the vertical velocities. However, the similar vertical velocity statistics do not indicate similar squall-line structures as shown in section 5.
c. Temperature, relative humidity, and precipitation comparisons
The Oklahoma Mesonet (hereinafter Mesonet) observations provide detailed surface data at 5-min intervals and allow for mesoscale comparisons between the observed and simulated squall-line evolution. To do so, the irregular station data from the Mesonet were interpolated onto a gridded domain with 1-km grid spacing covering most of Oklahoma. The interpolated temperature field at 9 m above ground level (AGL) at 0955 UTC is shown in Fig. 5a. The strong temperature gradient in the southern domain indicates the location of the squall-line leading edge. The low temperature area behind the line is the cold pool region. We adapt an Eulerian approach to compare the modeled and observed evolutions of the temperature, relative humidity, and precipitation fields within certain areas that are fixed in space.
For the observations, two rectangular areas of 120 km by 30 km that are aligned with the squall-line leading edge were chosen to calculate the time series of averaged values of the aforementioned fields in both the southern and the northern portions of the squall line (see Fig. 5a for locations of these two areas). These two areas illustrate the natural variability of the observed features. Both areas enclose at least three Mesonet sites so that the time series are not calculated based on interpolated data alone. The northern box is centered at Morris, Oklahoma, where the model sounding was taken. For the model results, the area with X between 320 and 350 km is used to calculate the corresponding time series for CNNB results (black box in Fig. 5b). Since the squall line moved faster in FSBM and UPNB than CNNB (the simulated average squall-line moving speeds after hour 2 are 16, 25, and 26 km h−1 for CNNB, FSBM, and UPNB, respectively while the actual squall line moved at about 31 km h−1 in the south and 20 km h−1 in the north), the area with X between 330 and 360 km is used in these two simulations (white box in Fig. 5b).
The 6-h time series of average low-level temperature, relative humidity, accumulated precipitation, and precipitation rate within the corresponding rectangular regions at 20-min intervals are illustrated in Fig. 6. A period of 6 h is needed to cover the entire evolution of the cold pool and precipitation features in the selected regions. For the southern and the northern Mesonet observations, data from 0800 to 1400 UTC and from 1100 to 1700 UTC are plotted, respectively. The observed temperature and relative humidity were measured at 9 and 2 m, while the model counterparts are from the lowest level (125-m height). Since we are primarily interested in the temperature drop and humidity increase during the passage of the squall line rather than the actual magnitudes of these quantities, for comparing the temperature time series (Fig. 6a), the northern observations were reduced by 1.7°C to match the model initial values. Similarly, the relative humidity of the northern observations was reduced by 3.7% (Fig. 6b). Because we initiated the model with the Morris sounding (center of the northern box), we anticipate that the temperature and moisture comparisons between model results and the northern box should be more representative. However, since many physical processes influencing the cold pool are not included in the current idealized model setup, accurate reproductions of the temperature and humidity features by the simulations are not necessarily expected.
As shown in Fig. 6a, temperature decreases of about 4°–5°C were observed as the squall line passed both the southern and the northern regions. CNNB simulates a weaker and slower temperature decrease of less than 3°C through 6 h. FSBM simulates a temperature trend similar to the northern observations but about 20 min slower. UPNB is able to reproduce the temperature evolution in the north during the first 3 h. Unlike the observations and the other two schemes, the UPNB-simulated temperature gradually increases after hour 3 and ends at similar values simulated by CNNB, which are about 1°C higher than the adjusted temperature observed in the north. Such a warming trend at low levels is a result of cold pool destruction by the strong descending rear inflow, which was not observed in the Mesonet data. More details about the adiabatic heating of the descending rear inflow will be discussed in section 5. Trends in the relative humidity are generally opposite of the temperature. The same trend from drier to moister conditions is seen for both the observations and simulations as the squall line passed. The slower moistening trend simulated by CNNB, relative to the other schemes, matches the northern observations well. FSBM simulates the most humid low-level atmosphere in the cold pool region corresponding to its strong cooling (Fig. 6a). UPNB captures the humidity trend in the southern area in the first 3 h. In response to the descending rear inflow in the simulations, the precipitation particles from the stratiform region completely evaporate before they reach the ground, causing a reduction of low-level relative humidity after hour 3.
The surface precipitation evolution in the CNNB simulation follows the northern observations in the first 3 h (Figs. 6c,d). As a result of the relatively slow speed of the CNNB-simulated squall line, the convective region stays in the study region longer than in the observations and with the other two schemes, and subsequently results in the highest accumulated precipitation among the simulations. FSBM simulates a similar precipitation pattern to the northern observations, while UPNB well captures the southern precipitation. FSBM is the only scheme that produces a stratiform precipitation rate close to the observations, in agreement with the reflectivity comparisons.
The Mesonet–model comparison together with the reflectivity and vertical velocity comparisons provide a consistent picture that the semi-idealized simulations using the three bin schemes capture many storm-scale features observed by a variety of instruments, but with notable differences in spatial and temporal patterns. However, significant intermodel differences occur. The cold pool dynamics and microphysics–dynamics interactions that lead to these scheme-specific features are analyzed and discussed in the next section.
5. Cold pool dynamics and dynamics–microphysics interactions
Interactions between the storm-generated cold pool and the low-level shear determine many storm-scale properties (RKW88; Weisman and Rotunno 2004). The idealized model setup provides a steady low-level wind shear and eliminates the impacts of radiation, land surface, and boundary layer processes on the cold pool evolution. Therefore, differences in the microphysics and the subsequent dynamical feedbacks are responsible for differences in the cold pool and squall-line evolution, simplifying the analysis and interpretation of simulation differences. The buoyancy perturbation B is defined following RKW88 as
where , , , and are potential temperature, vapor mixing ratio, total condensate mixing ratio (all hydrometeors are included), and gravity acceleration, respectively; and and by definition are the horizontally averaged potential temperature and vapor mixing ratio profiles, respectively. Here we use profiles from the initial sounding so that they are the same for all simulations and time independent. Term 1 is the temperature perturbation term, term 2 is the vapor perturbation term, and term 3 is the condensate loading term. In this study, the cold pool is defined as the region with buoyancy less than −0.01 m s−2 behind the convective core (see Figs. 10 and 12 for reference). The contour of B = −0.01 m s−2 adjacent to the convective region is defined as the cold pool leading edge. This threshold corresponds to a potential temperature perturbation of roughly −0.3 K, somewhat smaller than the threshold commonly used in literature (~−1 K) to define the cold pool (Bryan and Morrison 2012; Morrison et al. 2012, 2015). The location of the leading edge is insensitive to this threshold (within a few kilometers when this threshold varies from −0.01 to −0.07 m s−2).
a. General features of the cold pool
Vertical cross sections of the y-averaged B field overlaid with a few temperature contours and wind vectors are plotted in Fig. 7a for all simulations at hour 6. Only part of the domain with X between 200 and 500 km and Z below 12 km is shown. The corresponding horizontal cross sections at 2-km height indicating 3D flow patterns are shown in Fig. 7b. The cold pool is the area in green and blue colors in Fig. 7. Many common features occur in the three simulations. The negatively buoyant air behind the leading edge and the associated updraft in front of it are apparent in each simulation. A cold pool depth of around 4 km is associated with the height of the melting level. The updraft at the leading edge tilts from 2 km up, which indicates an unbalanced state between the vorticity associated with the cold pool and the ambient wind shear in which the cold pool dominates the shear (RKW88). As a result, the storm front inflow is almost horizontal above 4 km with limited net upward motion behind the leading edge. The dominance of the cold pool over the environmental shear also results in a descending rear inflow (Weisman 1992), which is apparent in each simulation.
Despite the similar buoyancy structures simulated by all of the schemes, the cold pool size, strength, moving speed, and the associated mesoscale flow patterns differ. Although Fig. 7b shows complicated flow patterns in the y direction, the overall dynamic and thermodynamic structure of the squall line can be described in 2D (x–z plane) reasonably well. Thus, the analysis is focused hereafter on y-averaged fields. Based on Fig. 7, CNNB generates the narrowest and the weakest cold pool, the slowest storm speed, and the most elevated and weakest rear inflow among the three schemes. FSBM simulates the widest and strongest cold pool, the strongest and deepest rear inflow, and a slightly slower storm moving speed than UPNB.
The cold pool intensity C is defined following Weisman and Rotunno (2004) as
where H is the cold pool depth, which is 4 km in this study, and is the horizontally averaged buoyancy within the cold pool. The value of C is a quantitative measure of the cold pool strength and the ratio of C to the wind difference between 0 and 4 km (the approximate height of the cold pool), dU, determines the relative balance between the vorticities associated with the cold pool and environmental shear (RKW88). Weisman (1992) also showed that the storm-generated rear inflow has a nonnegligible impact on the cold pool–shear balance. Therefore, the cold pool intensity should be adjusted for the presence of the rear inflow as
where C′ is the adjusted cold pool intensity and RI is the rear-inflow speed. Since the rear inflow generates a vorticity tendency that counteracts the vorticity tendency associated with the cold pool head, the cold pool intensity is reduced (Weisman 1992). In this study, C is calculated over the cold pool head, averaged over the region between the surface gust front and 30 km behind the gust front consistent with previous studies (Bryan and Morrison 2012; Lebo and Morrison 2014; Morrison et al. 2015). Calculations using other methods such as excluding positive buoyancy and varying the horizontal distance to be averaged did not show qualitative differences compared to the traditional approach used in this study. The term RI is represented by the maximum U at the back edge of the cold pool head.
Figure 8 illustrates the time series of C, RI, and C′ (Fig. 8a) and the ratios of C/dU and C′/dU in (Fig. 8b). These results show that FSBM simulated larger C than UPNB after hour 3, which agreed with Fig. 7. However, when the rear-inflow impact is considered, C′ is larger in UPNB for the first 5 h of simulations consistent with its faster squall-line moving speed than FSBM. CNNB generated the weakest cold pool and smallest C. The C/dU and C′/dU ratios are all greater than 1 indicating the dominance of the cold pool over the shear, resulting in tilted convective cores and descending rear inflow in all of the simulations. The rear inflow and its impacts on the storm dynamics and thermodynamics will be further discussed in section 5c.
b. Cold pool evolution and microphysics interaction
The evolution of the general storm structure depends on the microphysical process rates associated with phase change, the precipitation rate, and the precipitation efficiency (the latter defined as the ratio between average surface precipitation rate and the sum of the average condensation and deposition rates).
FSBM simulates the highest condensate production rate (condensation plus deposition), the highest evaporation plus sublimation rate and the lowest precipitation rate. The low densities and terminal velocities of ice-phase particles especially graupel in FSBM (see Fig. 1 for details) means that the solid particles experience a longer residence time before they completely melt or reach the ground and more ice mass is advected to the back of the storm. As the volume of the storm becomes larger, more ice mass falls into the stratiform region, which subsequently results in more deposition, sublimation, and evaporation than with the other schemes. The large condensate production rate in FSBM is consistent with somewhat stronger convective vertical velocities in FSBM (Figs. 4 and 11). It is possible that the balance between having the largest condensate production as well as the largest evaporation/sublimation rate in FSBM leads to less precipitation reaching the ground compared to other schemes. On the other hand, CNNB assumes high density for all ice-phase particles resulting in the highest fall speeds (Fig. 1). As a result, the lowest microphysical conversion rates and the highest precipitation rate are simulated by CNNB. The low condensation and deposition rates can also be attributed to CNNB having somewhat smaller vertical velocities (Figs. 4 and 11). UPNB assumes density and terminal velocity relationships between CNNB and FSBM (Fig. 1) and simulates intermediate rates. It appears that all schemes converge in terms of precipitation rate toward the end of the simulations.
The rate of mass change from freezing (melting) is comparable to vapor diffusion (Figs. 9a,b), with riming as the main process converting liquid water into ice. However, since the latent heat of fusion is about an order of magnitude smaller than the latent heats of condensation and deposition, the effect of freezing and melting on temperature is smaller. The net ice mass production rate above the melting layer (deposition + freezing − sublimation) is largest in UPNB, leading to the highest melting rate among the simulations. FSBM simulates a net ice production rate similar to UPNB, which in turn results in a higher melting rate than in CNNB (Fig. 9c). Precipitation rates and efficiencies from these simulations are consistent with the terminal velocity relationships, with the largest terminal velocity leading to the highest precipitation rate and the highest efficiency, and vice versa.
The cold pool strength is mostly determined by the temperature change because the temperature perturbation term dominates the other two terms in the cold pool buoyancy. The y-averaged temperature tendency due to microphysics is plotted in Fig. 10 for each scheme in the same subdomain as Fig. 7a at hour 6. Profiles of liquid water content (LWC), ice water content (IWC), vertical hydrometeor mass fluxes (downward flux as positive) in updrafts (w > 0 m s−1) and downdrafts (w ≤ 0 m s−1), and the average vertical velocity in updraft and downdraft regions are plotted in Fig. 11. For all three simulations, the profiles are calculated separately in the back of the storm (defined as the leftmost 30 km of the storm that experiences evaporation cooling; left panels in Fig. 11), behind the cold pool head (from 60 km behind to 30 km behind the cold pool leading edge; middle panels), and in the cold pool head (from 30 km behind to the cold pool edge; right panels). These three regions are indicated by purple dashed boxes in Figs. 10a1,b1,c1.
The strong and concentrated condensation heating in the convective cores and broad evaporative cooling in the stratiform region below the melting level are simulated by all schemes (Figs. 10a1,b1,c1). In the cold pool head region, FSBM simulates the smallest fall speeds and the strongest convective cores that lead to more upward transport of particles and more droplet and ice crystal activation (Fig. 11c). The evaporative cooling immediately upshear of the convective region is stronger in CNNB and UPNB than FSBM (Figs. 10a1,b1,c1) because of the higher ice particle fall speeds and thus larger hydrometeor fluxes into the melting layer (Fig. 11c2). However, because of the slower moving speed and smaller storm size, air parcels in rear-to-front flow at low levels within the cold pool spend less time experiencing evaporative cooling in CNNB than FSBM, which leads to its weaker cold pool head even though the evaporative cooling rate is stronger near the convective region in CNNB relative to the other simulations.
The larger IWC in the downdraft region above the melting layer in FSBM leads to stronger deposition and sublimation than in CNNB and UPNB (Figs. 10a2,b2,c2). Although the IWC and LWC are largest in FSBM above the melting layer (Fig. 11c1), it is speculated that its low riming efficiency makes the mass conversion rate from freezing (mostly riming) similar to the other two simulations (Figs. 10a3,b3,c3). Because of its larger graupel fall speeds, ice particles reach the ground in CNNB (Fig. 11c1) and it has the deepest melting layer (Fig. 10a3). All schemes simulate deeper melting layers in the cold pool head than in other areas due to the larger, denser, and faster graupel particles extending the melting process to the lower level in this region. Both CNNB and UPNB simulate larger hydrometeor mass fluxes in the downdraft region below the melting layer than FSBM, driving stronger evaporation in the cold pool head region (Figs. 10a1,b1,c1).
In the region behind the cold pool head, updrafts are much weaker than in the cold pool head region and little water exists above the melting layer (Fig. 11b). FSBM features the highest IWC due to its lowest fall speeds. Consequently, most IWC is transported away from the convective region (Fig. 11b1). As a result, FSBM simulates the strongest deposition, sublimation, and melting (Figs. 10a2–c2,a3–c3). The shallower and wider region of cooling farther upshear is the result of melting snow particles (mostly aggregates) that fall slower than graupel particles in the cold pool head region (Figs. 10a3–c3). The hydrometeor fluxes in the downdraft region below melting layer are in line with the evaporation rate patterns (Figs. 10a1–c1, 11b2).
In the back of the storm, mean upward motion is fairly weak, although the average downdraft around 2 km is stronger than in the other two regions for all of the schemes (Fig. 11a3). This region is where the rear inflow is generated (see section 5c). The highest IWC aloft is produced by FSBM, due to the lowest fall speed and active ice nucleation, leading to the strongest melting and sublimation rates among the simulations (Figs. 10b2,b3) similar to the other storm regions. This sublimation pattern generates the deepest descending rear inflow (Fig. 11a3). FSBM also has the largest mass flux and the strongest evaporation rate among the simulations in this region (Figs. 10c1 and 11a1,a2).
A trajectory analysis similar to Weisman (1992) is needed to quantify cold pool initiation and evolution. Limited by the 20-min outputs (the size of the output file is huge due to hundreds of 3D fields of all bin variables), we analyze the equivalent potential temperature θe, which is nearly conserved in precipitating moist adiabatic flow, to study how the cold pool evolves under the microphysics–dynamics interactions. The equivalent potential temperature (or moist static energy), as typically defined, is not conserved when ice processes are involved. However, due to sixfold difference between the latent heats of condensation and freezing, we do not expect the general features of the analysis to be significantly affected. Figure 12 shows three snapshots of the y-averaged θe field representing the initial (t = 20 min), developing (t = 100 min), and mature (t = 240 min) stages of the simulated squall lines. At the initial stage, warm and moist low-level high θe air is transported to high levels (Figs. 12a1–c1), with the convective updraft region splitting the cool and dry midlevel low θe air. The compensating downdraft region together with downdrafts induced by precipitation loading and evaporative cooling brings down low θe air behind the convective region toward the low levels.
At 100 min, the descending rear inflow develops and transports low θe air to the cold pool head region. Along the path, the air experiences adiabatic warming from descent and cooling from evaporation, sublimation, and melting (Fig. 10). The evolution of the cold pool air is a result of the balance between advection, adiabatic warming, and latent cooling. The returning surface flow also brings cold air toward the back of the domain. Low θe air in front of the line also contributes to cold pool development. As suggested by Fig. 12c2, midlevel environmental air penetrates the convective region and merges into the rear-inflow stream. This phenomenon was also found by RKW88. Convective updrafts are the most upright in CNNB due to its weak cold pool while they are strongly tilted in UPNB because of its strong cold pool (Fig. 8b), consistent with RKW theory. Interestingly, CNNB has the most upright convective drafts but weakest updraft velocities among all schemes (Figs. 9 and 11c3), which contradicts the consensus that more vertical updrafts should be stronger.2 However, in a strict sense the vorticity balance between the cold pool and environmental wind shear in RKW theory dictates the updraft geometry and tilt but not necessarily the strength (Weisman and Rotunno 2004; Bryan et al. 2006; Bryan and Rotunno 2014). After the squall line matures, the basic structure remains unchanged with updrafts tilting more upshear in FSBM and UPNB (Figs. 12b3,c3).
c. Rear-inflow generation and its effects on thermodynamics
As discussed previously, the storm-generated rear inflow is an important component of the squall-line structure that controls many features of the system. The rear inflow is associated with the vorticity tendency or the negative horizontal buoyancy gradient in the rear part of the stratiform region (Weisman 1992). Since the squall line is initiated around X = 300 km and there is no mean wind above 4 km, only the convective region moves forward and the back of the storm (tilted sublimation region) stays relatively steady within the purple box in Fig. 7a2. This area is defined as the rear-inflow generation region. The total vorticity tendency within this region is plotted as a function of time for each scheme in Fig. 13. Positive values indicate that positive vorticity tendency dominates, which implies a descending trend of the rear inflow [see details in Weisman (1992)]. The rear inflow strongly descends in all of the simulations before hour 2. Subsequently, it ascends for a brief period of 40 min and then stays relatively stable but with a weak descending tendency in CNNB. This ascending period leads to a rear inflow at higher altitude than in the other simulations. FSBM shows the strongest descending tendency overall and consequently simulates the lowest-altitude rear inflow. UPNB is between the other two simulations. The fastest rear-inflow speed is in FSBM and the slowest is in CNNB (Fig. 8a), also explained by the vorticity tendency in the rear-inflow generation region.
As a result of the descending rear inflow, strong evaporation is simulated by each scheme in the back of the storm (Figs. 10a1–c1). The low-level air in this region becomes warmer and dryer because the adiabatic warming of the descending air is not offset as much by evaporation as more of the hydrometeors have already evaporated at higher altitudes. Among the simulations, the least descending and slowest rear inflow in CNNB results in a small warm anomaly at X = 310 km and Z = 1 km (Fig. 7a1) and a relatively deep surface cold air layer caused by the deeper returning flow. At the other extreme, the most descending and fastest rear inflow in FSBM leads to a strong warm anomaly for X between 280 and 300 km and below 1 km. Such phenomena, called “heat bursts,” have been observed in real storms (e.g., Johnson et al. 1989), although they did not occur in this case based on the Mesonet data. The surface cold air layer behind the storm is blocked by the rear inflow in FSBM. UPNB simulates a surface cold air layer and warm anomalies between the extremes of CNNB and FSBM, consistent with its intermediate rear-inflow characteristics.
In attempt to show the sensitivity of simulation results to the “out of box” versions of different bin (spectral) microphysics schemes, the squall-line event on 20 May 2011, during the MC3E field campaign has been simulated by three bin microphysics schemes coupled into the WRF Model, version 3.4.1, in a semi-idealized setup. The 3D simulations driven by observed temperature and moisture profiles in the preconvection environment at 1200 UTC in Morris, Oklahoma, show that each scheme produced a squall line with features broadly consistent with the observed storm characteristics (e.g., reflectivity pattern, vertical velocity statistics, and evolution of low-level air temperature, relative humidity, and precipitation). However, substantial differences in the details of the system are evident among the simulations. For instance, the reflectivity patterns are remarkably different (cf. Fig. 3), cold pool extent and strength vary dramatically (cf. Fig. 7), and the atmospheric thermodynamic and flow structure of the system deviates among the simulations (cf. Fig. 12). Because of the idealized setup of the experiments, all differences are attributable to the different representations of microphysical processes and properties, and the resulting interactions between the microphysics, cold pool, and dynamics.
CNNB assumes high densities for ice-phase particles resulting in the highest fall speeds among all of the schemes tested, which led to rapid fall out and less time for precipitation particles to melt, sublimate, and evaporate. As a result, the highest precipitation rate and the weakest cold pool head are simulated by CNNB. Weak storm-relative inflow associated with the weak cold pool and the large hydrometeor fall speeds limited the upshear transport of ice particles aloft. The small amount of ice aloft resulted in weak sublimation in the back of the storm and hence weak descending rear inflow in CNNB. In contrast, FSBM produced the highest condensate production and evaporation and the lowest precipitation efficiency primarily due to its low terminal velocities associated with the density-size assumptions of ice particles especially graupel. As a result of its strong evaporation, FSBM produced the strongest cold pool. Large sublimation rates in the back of the storm generated a strong descending rear inflow, which led to strong heating near the ground. The vorticity tendency associated with the strong descending rear inflow counteracted the vorticity tendency associated with the cold pool head, reducing the storm moving speed (Weisman 1992). UPNB assumes ice-phase particle densities and terminal velocities between CNNB and FSBM and simulated condensate production, evaporation, precipitation efficiency, cold pool, and rear-inflow strength between the other two simulations. The balance between the rear inflow and the cold pool strength leads to the fastest squall-line moving speed in UPNB. The terminal velocity relationship has a great effect on simulated precipitation. Indeed, the collision efficiency and diffusion rate depend on such a relationship through the relative fall velocity and ventilation process. However, the different precipitation rates simulated by these schemes are not just the results of different fall speed relationships. Other mechanisms contribute to such differences as well. For instance, the instant melting assumed by FSBM leads to stronger evaporation initially and moister lower-level atmosphere, which suppresses further evaporation and enhances precipitation rate at the later stage.
The bin microphysics intercomparison study presented here used the same dynamical framework with physical processes besides microphysics turned off and thus the differences in model results are solely due to differences in the microphysics schemes and their interaction with the dynamics. This study shows that bin microphysics schemes, designed to be conceptually more realistic and thus arguably more accurate than bulk microphysics schemes, simulate a wide spread of microphysical, thermodynamic, and dynamic characteristics of a squall line even in this constrained modeling framework. Users of these bin schemes should be clear about what detailed hydrometeor properties are assumed within the schemes to effectively apply them to their specific research topics. For instance, the low density and small terminal velocity of the graupel particles in the default version of FSBM may not correctly represent the hydrometeor properties in the convective cores of a squall line. When a version of FSBM used by Fan et al. (2017) that employs hail particles for high-density ice particles was used to simulate the same case, the simulated squall line showed a much stronger reflectivity field in the core region, in a better agreement with observations. This result reinforces our conclusion that the simulated squall-line structure is sensitive to the fundamental assumptions of the ice particle properties.
Differences in the squall-line characteristics simulated by the three bin schemes in this study are qualitatively similar to the differences using various bulk microphysics schemes in Morrison et al. (2015), with the caveat that the case in this study is different from the squall-line case in Morrison et al. (2015). Many of the differences produced by the bin and bulk schemes are attributed to similar causes, such as narrower (wider) convective regions with higher (lower) reflectivity simulated by schemes assuming ice particles with high (low) terminal velocity. Given this uncertainty in the treatment of microphysical processes and ice particle properties in particular, it is recommended to constrain bin simulations with detailed observations for the use of developing and testing bulk schemes for deep convective cases. Future work may focus on improving the representation of ice particle properties in bin schemes to reduce this uncertainty and using similar assumptions for all schemes to isolate the impact of physics from numerics. One possible avenue for such work would be to move away from the paradigm of ice categories with fixed particle properties and instead use categories with predicted, freely evolving properties (e.g., Hashino and Tripoli 2007; Harrington et al. 2013; Morrison and Milbrandt 2015).
This research was partially supported by the U.S. Department of Energy (DOE) Atmospheric System Research (ASR) Program (DE-SC0008648 and DE-SC0016476) and was partially supported by the National Science Foundation, through funding to the National Center for Atmospheric Research (NCAR). The National Center for Atmospheric Research is sponsored by the National Science Foundation. J. Fan’s contribution was supported by the DOE ASR Program. The Pacific Northwest National Laboratory is operated for the DOE by Battelle Memorial Institute under Contract DE-AC06-76RLO1830. Z. Lebo was partially supported by the National Oceanic and Atmospheric Administration’s Climate Goal funding and University of Wyoming startup funding. W. Wu and G. McFarquhar appreciate the support of the office of Biological and Environmental Research (BER) of the DOE (DE-SC0014065) and the support from the NCAR Advanced Study Program (ASP) faculty visiting program. X. Chu was supported by the NCAR ASP graduate student visiting program. The contribution to this research by I. Geresdi was supported by the Hungarian Scientific Research Fund (Development and application of novel numerical model to investigate the precipitation formation in mixed phase clouds). X. Lou was supported by the National Science Foundation of China (41275148). The authors appreciate the help from Drs. Xiquan Dong, Pavlos Kollias on radar data analysis, the scripts from Dr. Noemi Sarkadi for detailed radar reflectivity calculation, and the helpful discussions with Drs. George Bryan and Morris Weisman.
In such simulations, the model domain and physics are set up in an idealized way but the initial thermodynamic profile is based on an observed sounding.