A new two-moment cloud microphysics scheme predicting the mixing ratios and number concentrations of five species (i.e., cloud droplets, cloud ice, snow, rain, and graupel) has been implemented into the Weather Research and Forecasting model (WRF). This scheme is used to investigate the formation and evolution of trailing stratiform precipitation in an idealized two-dimensional squall line. Results are compared to those using a one-moment version of the scheme that predicts only the mixing ratios of the species, and diagnoses the number concentrations from the specified size distribution intercept parameter and predicted mixing ratio. The overall structure of the storm is similar using either the one- or two-moment schemes, although there are notable differences. The two-moment (2-M) scheme produces a widespread region of trailing stratiform precipitation within several hours of the storm formation. In contrast, there is negligible trailing stratiform precipitation using the one-moment (1-M) scheme. The primary reason for this difference are reduced rain evaporation rates in 2-M compared to 1-M in the trailing stratiform region, leading directly to greater rain mixing ratios and surface rainfall rates. Second, increased rain evaporation rates in 2-M compared to 1-M in the convective region at midlevels result in weaker convective updraft cells and increased midlevel detrainment and flux of positively buoyant air from the convective into the stratiform region. This flux is in turn associated with a stronger mesoscale updraft in the stratiform region and enhanced ice growth rates. The reduced (increased) rates of rain evaporation in the stratiform (convective) regions in 2-M are associated with differences in the predicted rain size distribution intercept parameter (which was specified as a constant in 1-M) between the two regions. This variability is consistent with surface disdrometer measurements in previous studies that show a rapid decrease of the rain intercept parameter during the transition from convective to stratiform rainfall.
Squall lines with trailing stratiform precipitation are common in both midlatitude and tropical environments, and have been extensively studied by numerous researchers (e.g., Fujita 1955; Zipser 1969; Ogura and Chen 1977; LeMone et al. 1984; Houze 1989; Biggerstaff and Houze 1991, 1993). These studies have suggested several common morphological features described by the conceptual model of Biggerstaff and Houze (1991, see their Fig. 18). These features include an upshear-tilted, multicellular convective region with heavy precipitation and active updraft cell generation along the gust front and a low-level radar reflectivity minimum between the convective and stratiform regions, followed by a region of moderate precipitation in the trailing stratiform region.
The kinematic structure is described in detail by Biggerstaff and Houze (1991) and references therein. The trailing stratiform region is associated with a mesoscale updraft at mid- and upper levels, and mesoscale downdraft on the order of tens of centimeters per second. Modeling studies (e.g., Fovell and Ogura 1988; Tao et al. 1993) and thermodynamic retrievals (e.g., Hauser et al. 1988) indicate that air at mid- and upper levels of the trailing stratiform region is positively buoyant, at least in part due to the advection of positively buoyant air from the convective region in the front-to-rear relative flow at mid- and upper levels. Several studies (e.g., Smull and Houze 1985; Biggerstaff and Houze 1991; Rutledge and Houze 1987) have suggested the important role of advection of hydrometeors from the convective region in generating trailing stratiform precipitation. However, additional growth of hydrometeors within the stratiform regions may also be important. Based on the results of a diagnostic modeling study, Rutledge and Houze (1987) found that only one-fourth of the stratiform rainfall would have reached the surface without additional particle growth generated by mesoscale ascent in the stratiform region.
The mesoscale downdraft appears to be associated with evaporation of precipitation falling into unsaturated air beneath the stratiform cloud (e.g., Brown 1979), which may be enhanced by melting (Leary and Houze 1979). Rutledge et al. (1988) found that the mesoscale downdraft was strongest just beneath the melting layer of the stratiform region in the area of heaviest precipitation. In a model sensitivity study, Tao et al. (1995) also found that melting (in both the stratiform and convective regions) had an impact on convective characteristics, producing upshear-tilted multicellular rather than unicellular structure. The unicellular structure was associated with more intense updrafts and reduced rearward horizontal mass flux at mid- and upper levels. Ferrier et al. (1996) also found that a transition from upshear-tilted multicellular to unicellular convection decreased the rearward flux of condensate and hence reduced the intensity of the stratiform rain region.
Many squall lines with trailing stratiform precipitation exhibit a distinct minimum in low-level radar reflectivity between the stratiform and convective regions. Several explanations for this transition zone have been put forth. Observational studies have suggested enhanced subsidence in the transition zone associated with dynamically and/or microphysically driven downdrafts (i.e., by evaporative cooling and precipitation loading; e.g., Houze and Rappaport 1984; Smull and Houze 1985; Srivastava et al. 1986). The diagnostic modeling study of Rutledge and Houze (1987) suggested that the location of the heaviest stratiform precipitation could be explained by the combination of upper-level storm-relative flow and slow hydrometeor fall speeds for particles originating from the tops of the convective cells. The importance of fall speed sorting on the distribution of precipitation was also suggested by the modeling study of Fovell and Ogura (1988). Biggerstaff and Houze (1993) argued that the transition zone was the result of a lack of crystal aggregation, in contrast to the stratiform region with a steady supply of small crystals available for aggregation. The lack of small crystals in the transition zone was presumably the result of sublimation due to midlevel subsidence embedded within the overall front-to-rear midlevel ascent.
Numerous studies have attempted to simulate squall lines with an enhanced region of trailing stratiform precipitation using cloud-system-resolving models (e.g., Fovell and Ogura 1988; Lafore and Moncrieff 1989; McCumber et al. 1991; Ferrier et al. 1995). Fovell and Ogura (1988) and McCumber et al. (1991) found that inclusion of ice microphysics enhanced stratiform precipitation relative to liquid-only schemes, because of an increased rearward flux of hydrometeors in the presence of snow/graupel (which tends to have a lower mean terminal fall speed than rain). However, the region of stratiform precipitation was still narrower and less intense than observed, and Fovell and Ogura (1988) further noted the lack of a clear transition zone of low radar reflectivity in their simulation as was observed. Subsequent studies (Sui et al. 1998; Lang et al. 2003) focused on the ratio of stratiform and convective precipitation; these studies also suggested a bias toward heavy (convective) rain rates, although Lang et al. (2003) note that their results were sensitive to the method used to partition stratiform and convective rain. In general, models have not been very successful at replicating commonly observed squall-line features, the transition zone, and the trailing stratiform region in particular.
Microphysics parameterizations in most cloud models are bulk schemes; that is, they assume an underlying shape for the hydrometeor size distribution, and predict one or more bulk quantities of the distribution. Detailed bin microphysics schemes that explicitly predict evolution of the size distribution are computationally demanding and therefore not feasible for most applications. An improvement in bulk microphysics schemes has been the prediction of two moments of the hydrometeor size spectra rather than just one (i.e., the mass mixing ratio). Several two-moment schemes have been developed and utilized in a variety of applications (e.g., Koenig and Murray 1976; Ferrier 1994; Meyers et al. 1997; Khairoutdinov and Kogan 2000; Seifert and Beheng 2001; Milbrandt and Yau 2005; Morrison et al. (2005, hereafter MCK). These schemes predict the number concentrations and mixing ratios of the hydrometeor species, which increases the degrees of freedom and potentially improves representation of the particle size distributions.
In the current study, a new two-moment microphysics scheme is implemented into the Weather Research and Forecasting model (WRF; Skamarock et al. 2007). The purpose of this paper is to document the new scheme, and to describe the impact of this scheme on the evolution of an idealized two-dimensional squall line compared to a one-moment version of the same scheme, with a focus on the evolution of the trailing stratiform region. While several studies have examined sensitivity of midlatitude and tropical squall lines to different microphysics schemes in terms of broad categories such as number or type of ice classes (no ice versus multiple-class ice schemes; e.g., Fovell and Ogura 1988; McCumber et al. 1991; Ferrier et al. 1995; Liu and Moncrieff 2007), to our knowledge no studies have focused specifically on the impact of predicting one versus two moments of the particle size distributions using different versions of the same scheme in terms of squall-line organization and structure. This issue may be especially important because two-moment schemes are beginning to be more widely used in models. Here we use a 2D framework similar to many previous studies (e.g., Fovell and Ogura 1988; Lafore and Moncrieff 1989; McCumber et al. 1991; Ferrier et al. 1995; Tao et al. 1995). The justification is that variability perpendicular to the squall line tends to be much greater than variability parallel to the line, at least over larger scales.
2. Description of the microphysics schemes
a. Two-moment scheme
The new two-moment microphysics scheme implemented into WRF predicts the mass mixing ratios and number concentrations of five hydrometeor species: cloud droplets, cloud ice, snow, rain, and graupel. This scheme is based on the parameterization of MCK, and subsequently implemented into the fifth-generation Pennsylvania State University–National Center for Atmospheric Research (PSU–NCAR) Mesoscale Model (MM5; Morrison and Pinto 2005, 2006; Morrison et al. 2008). The most significant difference relative to the earlier scheme of MCK is the addition of prognostic variables for graupel mixing ratio and number concentration. Key aspects of the scheme as well as other differences with respect to MCK are described below.
The cloud and precipitation particle size distributions are represented by gamma functions:
where N0, λ, and μ are the intercept, slope, and shape parameters of the size distribution, respectively, and D is the particle diameter. The parameters N0 and λ are derived from the predicted number concentration N and mixing ratio q, and the specified μ for each species:
where Γ is the Euler gamma function and the parameters c and d are given by the assumed power-law mass–diameter (m–D) relationship of the hydrometeors for each species, where m = cDd. Here, all particles are assumed to be spheres for simplicity, with a bulk particle density for the various ice species given by Reisner et al. (1998). For the precipitation species (i.e., rain, snow, and graupel), as well as cloud ice, μ = 0. Thus, the size distributions for these species are exponential functions (Marshall–Palmer distributions). For cloud droplets, μ is a function of the predicted droplet number concentration following the observations of Martin et al. (1994). In MCK, μ for cloud droplets was derived from the theoretical formulation of Khvorostyanov and Curry (1999).
Equations for the time tendencies of cloud droplets, cloud ice, snow, and rain mixing ratios and number concentrations are given by MCK [see their Eqs. (1) and (2)]. The additional time tendency equations for graupel mixing ratio and number concentration are given by Reisner et al. (1998), except that minimum mixing ratios of 0.1 g kg−1 for rain and snow are required to produce graupel from collisions between rain and snow; minimum mixing ratios of 0.1 and 0.5 g kg−1 for snow and cloud water are required to produce graupel from collisions between snow and cloud droplets; and the minimum mixing ratio of 0.1 g kg−1 for rain is required to produce graupel from collisions between rain and cloud ice, following Rutledge and Hobbs (1984). The threshold mixing ratios for conversion to graupel are fairly arbitrary; sensitivity to these thresholds was examined by Morrison and Grabowski (2008). Further investigation of the impact of the graupel microphysics is beyond the scope of this study.
The change in snow and graupel number concentration due to melting and sublimation (not included in Reisner et al. 1998) is calculated by assuming that the relative change in number concentration is the same as the change in mixing ratio due to these processes following Ferrier (1994). Following MCK and Ferrier (1994), we assume that during rain evaporation the relative decrease of number concentration is the same as mixing ratio, which implies that the mean size does not change (the decrease of number concentration during sublimation of snow and graupel is similarly treated). The increase in rain number concentration due to melting is equal to the decrease of graupel and snow number concentration. Note that these processes are not treated in one-moment schemes since they do not predict particle number concentrations.
Radar reflectivity, Ze, is calculated from integration of the size distributions for each species following Smith (1984). For simplicity, and since we are primarily interested in comparisons to widely available radar data, only Rayleigh scattering is considered whereas Mie scattering is ignored. This assumption is justified for the relatively large wavelength used in this study (10 cm). For frozen species, a prefactor is used to compensate for the fact that the dielectric factor is with respect to water not ice. The special case of partially melted snow and graupel utilizes code of Blahak (Blahak 2007), which allows for different ice lattice and water coating assumptions. This produces a radar bright band that appears physically reasonable and improves upon the assumption of no meltwater.
b. One-moment scheme
The one-moment scheme is exactly the same as the two-moment scheme described above, except that the number concentrations of the precipitation species are diagnosed, not predicted. Thus, N0 is specified for rain, snow, and graupel, and N and λ are derived from the predicted q and specified N0 and μ for each species by rearranging terms in (2) and (3). Following many one-moment schemes (e.g., Lin et al. 1983; Rutledge and Hobbs 1984; Dudhia 1989; Grabowski 1998), N0 is constant for a given species. Values are taken from existing schemes described in the literature and shown in Table 1. Note that while most one-moment schemes use constant values of N0, it is not an intrinsic feature of these schemes; some of them allow N0 to vary for rain, snow, and/or graupel as a function of the predicted cloud or thermodynamic variables (e.g., temperature or mixing ratio of the species). For example, Reisner et al. (1998) and Hong et al. (2004; i.e, WSM3, WSM5, and WSM6 schemes in WRF) allow variable N0 for snow, while Thompson et al. (2004, 2008) allow variable N0 for all precipitation species (i.e., snow, graupel, and rain).
3. Experimental design
The Advanced Research WRF model, version 2.2 (Skamarock et al. 2007) is the fully compressible, nonhydrostatic, two-dimensional cloud model used for the simulations. Governing equations are solved using a time-split integration with third-order Runge–Kutta scheme and substeps for the acoustic and gravity wave modes. Prognostic variables include the 2D velocity components, perturbation potential temperature, perturbation geopotential, perturbation surface pressure of dry air, water vapor mixing ratio, and the various cloud microphysics variables.
The setup for this case follows from the standard WRF 2D idealized squall-line case available as part of the WRF modeling system. Lateral boundary conditions are open. Horizontal and vertical turbulent diffusion are calculated using a 1.5-order turbulent kinetic energy (TKE) scheme (Skamarock et al. 2007). The vertical domain extends from the surface to 20 km. A Rayleigh damper with damping coefficient (inverse damping time scale) of 0.003 s−1 is used in the top 5 km to damp spurious waves in the stratosphere. The upper and lower boundaries are free slip with zero vertical velocity. Surface fluxes are zero and radiative transfer is neglected. Previous studies have shown some impact of radiative transfer on squall-line development (Tao et al. 1993, 1996), but these impacts are neglected here for simplicity.
The model is initialized with the environmental temperature and moisture profiles used by Weisman and Klemp (1982, 1984) and others (see Fig. 1 in Weisman and Klemp 1982), which are broadly typical of the environment of midlatitude, continental squall lines although they are relatively moist at midlevels. The initial convective available potential energy (CAPE) is 2200 J kg−1. The melting level is located at a height of 4 km. The horizontal wind profile has a shear of 0.0048 s−1 in the lowest 2.5 km and zero shear above. The mean wind above 2.5 km is zero, which helps to keep the squall line near the center of the domain. Convection is initially triggered by adding a thermal with maximum perturbation in potential temperature of 3 K centered at a height of 1.5 km and varying as the cosine squared to the 0-K perturbation at its edge. The thermal has a horizontal radius of 4 km and a vertical radius of 1.5 km. We also ran simulations with convection triggered by adding a cold pool in a portion of the domain (−5 K perturbation potential temperature near surface decreasing linearly to zero at 3 km); results are qualitatively similar to those using a thermal and are therefore not shown here. The model uses 1000- and 250-m horizontal and vertical grid spacing, respectively, over a domain that is 600 km in the horizontal. The model time step is 5 s, with a simulation period of 7 h. Results show little sensitivity to increasing the domain size. Reducing the horizontal grid spacing to 250 m changes details of the simulation, but the overall structure is similar. Note that Bryan et al. (2003) suggested that grid spacing of at least o(100 m) may be needed for traditional subgrid turbulence closures in large-eddy models when simulating deep convection in three dimensions.
Results using both one-moment (1-M) and two-moment (2-M) schemes are analyzed in detail. In addition, several sensitivity tests are described in the next two sections. The results shown here are based on model output archived every 10 min of simulation time.
a. The 2-M simulation
Moist deep convection is initiated within the first 15 min of the simulations and, over the next several hours, it organizes into features that are characteristic of many observed squall lines. The evolution of the storm in 2-M is shown by time series of maximum and minimum vertical velocity, minimum perturbation potential temperature, θ′, and domain-average rain rate (Fig. 1). In addition, plots of radar reflectivity, Z, and 2D wind vectors are shown at 4 and 6 h (Figs. 2 –3). Updrafts and downdrafts are stronger and more transient in terms of peak intensity during the first 3 h of the simulation, before reaching a mature quasi-equilibrium phase. Domain-average rain rate increases and θ′ decreases during the initial phase before reaching quasi-equilibrium after approximately 3 h of integration.
By 4 h, an accompanying trailing stratiform region has appeared (Fig. 2) after transitioning from a more symmetric initial phase. Well-developed front-to-rear flow at mid- and upper levels, as well as rear-to-front flow at mid- and lower levels, is consistent with observed squall lines (e.g., Biggerstaff and Houze 1991, their Fig. 18). At lower levels, the rear-to-front flow extends to the surface. The updraft cells exhibit weak upshear tilt. These cells move away from the leading edge in the front-to-rear flow at mid- and upper levels and subsequently weaken, thereby allowing rapid fallout of precipitation and intensification of low-level downdrafts beneath the cell. New cells periodically arise along the gust front. Once this mature stage is fully developed within about 3–4 h, the average forward propagation speed (grid relative) of the gust front is about 6 m s−1 (keeping in mind that the ambient wind above the low-level shear is zero). Near-surface temperatures cool by as much as ∼10 K behind the gust front.
By 6 h, the trailing stratiform region is fairly well defined, with the region of stratiform precipitation at the surface extending over a region 130 km wide (Fig. 3). This region is associated with mesoscale updrafts and downdrafts above and below the melting layer (Fig. 4). The mesoscale updraft and downdraft are similar in magnitude (∼30–60 cm s−1). The strongest part of the mesoscale downdraft is shifted toward the convective region relative to the mesoscale updraft, with peak intensity just below the melting layer, consistent with observations (e.g., Rutledge et al. 1988; Biggerstaff and Houze 1991, 1993). However, a potential weakness of the simulation is the lack of a clearly defined transition region of lighter precipitation between the convective and stratiform regions, which is often observed (e.g., Biggerstaff and Houze 1993).
b. Comparison of the 1-M and 2-M simulations
1) Storm morphology
The overall storm morphology is similar between 1-M and 2-M (i.e., a propagating squall line with a sharp leading edge of convective precipitation followed by a trailing stratiform region). The 1-M and 2-M results are quite similar over the first 2–3 h during the more symmetric initial phase of the storm. However, by 6 h, the simulations have noticeably diverged in some aspects as shown by a comparison of Figs. 3 and 5. The most significant difference is the much narrower and weaker region of trailing stratiform precipitation in 1-M compared to 2-M. Radar reflectivities in 1-M are about 5–20 dBZ lower than 2-M both above and below the melting layer in the stratiform region. Above the melting layer, the region of high reflectivity values (>30 dBZ) extends about 40 km farther behind the storm’s leading edge in 2-M. Differences are even more noticeable below the melting layer, with very little rain reaching the surface in 1-M outside of the core convective region. The rapid decrease of radar reflectivity between the melting layer and the surface in the stratiform region in 1-M is indicative of large rain evaporation rates, which is described in more detail later. In contrast, radar reflectivity is nearly constant with height in the stratiform region in 2-M below the melting layer due to lower evaporation rates, which is more consistent with observed reflectivity profiles (e.g., Biggerstaff and Houze 1993, their Fig. 9). The profiles of radar reflectivity in 2-M are also impacted by raindrop size sorting (an effect that is absent in 1-M). This mechanism produces relatively large reflectivities near the surface associated with large mean drop size compared with smaller drop size higher in the cloud.
2) Precipitation and the cold pool
Differences in the trailing stratiform regions between 1-M and 2-M are further revealed by comparing surface rainfall rates. Rain rates in the stratiform region are much larger in 2-M, while they tend to be slightly smaller in the convective region (Fig. 6). However, the domain-average precipitation rate is fairly similar between the runs (see Fig. 1d). Interestingly, superimposed on the quasi-equilibrium behavior in 1-M (i.e., the periodic generation of new cells along the gust front) is a period of discrete propagation due to initiation of convection well ahead (by a few tens of kilometers) of the gust front at about 3 h (seen by the sharp jog of the leading edge of precipitation at X ≈ 370 km in 1-M in Fig. 6). This discrete propagation is insensitive to the domain size and is caused by midlevel (approximately 4–8 km) vertical velocity perturbations that are likely associated with vertically trapped gravity waves. See Fovell et al. (2006) for a detailed discussion of this mechanism.
To illustrate further differences between 1-M and 2-M, we examine the average rain mixing ratios and evaporation rates between 6 and 7 h, after the simulations have attained quasi-equilibrium (Figs. 7 and 8). Results in general are insensitive to the particular time period analyzed (except during the episode of discrete propagation in 1-M). As shown in Fig. 8, rain evaporation rates are several times larger in 1-M, especially just below the melting layer and in the forward portion of the stratiform region near the convective region (at X ≈ −40 to −150 km in Fig. 8). This enhanced evaporation leads to rain mixing ratios that are smaller in 1-M than 2-M near the surface in the stratiform region (see Fig. 7). Rain mixing ratios and evaporation rates are larger in 2-M than 1-M at midlevels in the convective region. This appears to have important consequences for the convective drafts as described later in this section.
Differences in the rain evaporation rate between 1-M and 2-M are associated with differences in the rain size distribution parameters (the rain intercept parameter N0r and the rain slope parameter λr) between the runs. The rain evaporation rate is given by (similar to Rutledge and Hobbs 1983; Reisner et al. 1998)
where Dυ is the diffusivity of water vapor in air, S is the liquid water saturation ratio, A′ and B′ are thermodynamic parameters related to the release of latent heat, f1 and f2 are ventilation parameters, ar and br are fall speed parameters for rain (where fall speed as a function of particle diameter D is given by arDbr), μa is the dynamic viscosity of air, Sc is the Schmidt number, and Γ is the Euler gamma function. In 2-M, N0r and λr are derived from the predicted rain mixing ratio and number concentration. In 1-M, N0r is specified here as a constant, and λr is derived from the specified N0r and the predicted mixing ratio (see section 2). For a given rain mixing ratio, evaporation rate depends on N0r only (in terms of the dependence on the size distribution parameters) since λr can be expressed in terms of N0r by combining and rearranging (2) and (3). Hence, we express differences in rain evaporation rate between 1-M and 2-M mostly in terms of differences in N0r. However, it should be kept in mind that differences in the rain mixing ratio between 1-M and 2-M also lead to differences in evaporation rate; thus, for a given value of N0r, larger rain mixing ratio will lead to a greater evaporation rate.
The predicted N0r in 2-M, averaged between 6 and 7 h, is shown in Fig. 9. A clear trend is evident, with the largest values of N0r occurring in the convective region, and steadily decreasing through the trailing stratiform region. This trend is due to significant droplet collision–coalescence in the convective region associated with large amounts of cloud water, which produces a large rain number concentration and correspondingly large N0r. In contrast, there is little droplet collision–coalescence in the stratiform region (as rain is produced almost entirely by melting of snow and graupel), resulting in relatively low rain number concentration and smaller N0r. The N0r and λr also tend to increase with height in both regions, which mostly reflects size sorting (i.e., separation of number concentration and mixing ratio due to the different number- and mass-weighted sedimentation velocities in 2-M). Note that two-moment schemes with exponential size distributions may produce excessive size sorting (Wacker and Seifert 2001). However, here the implicit parameterization of rain drop breakup [by limiting λr to a minimum of 20 cm−1 following Hodson (1986), Srivastava (1987), and others] largely controls the distribution of λr near the surface in 2-M, which prevents unrealistically large mean particle size due to excessive size sorting. Diagnostically varying the shape parameter μ (Milbrandt and Yau 2005; Seifert 2008) to allow for nonexponential size distributions may also curtail excessive size sorting.
In most of the trailing stratiform region, the predicted N0r is between 105 and 107 m−4 (except for larger values near the surface that are associated with droplet collision–coalescence in a thin, low-level cloud that forms in the cold pool). The mean value is 2 × 106 m−4 (for all points X < −100 km in Fig. 9 with rain mixing ratio > 0.01 g kg−1, excluding the large near-surface values), compared to the specified value of 107 m−4 in 1-M. Thus, the predicted N0r in 2-M is associated with an evaporation rate that is approximately a factor of 2 smaller than in 1-M in the stratiform region, for given environmental conditions and rain mixing ratios. The smaller values of N0r in 2-M than 1-M in the stratiform region are also consistent with smaller λr, which directly contributes to larger radar reflectivity for a given rain mixing ratio because of the larger mean rain drop size associated with the smaller λr.
In the convective region the N0r predicted in 2-M is generally larger than 107 m−4 (approaching 109 m−4 near the melting layer), except near the surface. The mean value in the convective region (for all points between −50 < X < 0 km in Fig. 9 with rain mixing ratio >0.01 g kg−1) is 2 × 108 m−4, compared to the specified value of 107 m−4 in 1-M. In contrast to the stratiform region, the larger values of predicted N0r in 2-M in the convective region at midlevels are associated with the larger rain evaporation rate there compared to 1-M.
The minimum perturbation potential temperature does not greatly differ between 1-M and 2-M (generally <1 K) once the quasi-equilibrium phase is reached after about 3 h of integration (see Fig. 1c). However, the increased evaporation rate in the stratiform region leads to a broader and generally colder cold pool in 1-M than 2-M, by up to about 3 K at low levels in the region between X ≈ −150 to −40 km (Fig. 10). This is consistent with a slightly increased propagation speed of the gust front in 1-M than 2-M (by an average of 1.8 m s−1), excluding the period of discrete propagation mentioned previously. Note that near the head of the cold pool (just behind the gust front), perturbation potential temperature is similar or even slightly colder in 2-M than 1-M.
Figure 11 shows that total ice mixing ratios (i.e., cloud ice, snow, and graupel) in the stratiform region are similar between 1-M and 2-M. However, earlier in the mature phase of the simulation (i.e., between 4 and 6 h), the stratiform region in 2-M has much greater amounts of ice condensate than 1-M, which also contributes to greater amounts of rain below the melting layer in the stratiform region in 2-M as the ice precipitation falls and melts. Conversely, the amount of ice condensate is much smaller in 2-M than 1-M in the convective region.
3) Mesoscale dynamics
Analysis of the mean vertical motion in the stratiform region is complicated by its small magnitude and large variability, presumably due to large-amplitude gravity waves. Nonetheless, useful patterns emerge with appropriate temporal and spatial averaging. A comparison of vertical velocities averaged between 6 and 7 h and over 20-km horizontal segments reveals that 2-M has a weaker downdraft (by a maximum of ∼30 cm s−1) at low levels of the stratiform region (Fig. 12) consistent with the reduced evaporative cooling described previously. In contrast, the mesoscale updraft at midlevels (∼4–7 km) is a few tens of centimeters per second stronger in 2-M than 1-M.
To further investigate differences between 1-M and 2-M in terms of the mesoscale updraft, we examine the latent heating rate and rearward horizontal fluxes of condensate and buoyancy1 averaged between 6 and 7 h. Air in the stratiform region at midlevels (∼4–6 km) in 2-M is more positively buoyant than it is in 1-M (Fig. 13), leading to a greater buoyancy gradient across the stratiform region. Two factors contribute to this difference in buoyancy. First, the latent heating rate due to vapor depositional growth of ice, and hence contribution to buoyancy via temperature perturbation, is somewhat greater in 2-M than 1-M in the stratiform region (Fig. 14). Second, there is a much larger flux of buoyancy at midlevels from the convective to the stratiform region in 2-M than 1-M (Fig. 15a). Greater updraft strength and penetrative depth of the convective cells in 1-M, due mostly to the different treatment of rain microphysics [detailed in section 4b(4)], leads to the reduced detrainment of positively buoyant air at midlevels (but increased detrainment at upper levels, explaining the larger upper-level rearward fluxes of buoyancy in 1-M as seen in Fig. 15a). A similar result was described by Tao et al. (1995) and Ferrier et al. (1996), who found that increased convective updraft strength reduced the detrainment and horizontal transport of buoyancy and momentum at midlevels.
In contrast to the buoyancy, the rearward flux of condensate from the convective region tends to be much smaller in 2-M than 1-M (Fig. 15b). This decreased flux is a result of the much smaller amounts of ice condensate in the convective region in 2-M than 1-M (see Fig. 11). Thus, the larger ice growth rates in the stratiform region in 2-M, associated with the stronger mesoscale updraft, compensate for the much smaller flux of condensate from the convective region, and result in similar amounts of ice in the stratiform region compared to 1-M as shown in Fig. 11 (and larger amounts in 2-M than 1-M during the early part of the mature stage of the squall line).
4) Convective drafts
Maximum updraft and downdraft velocities tend to be greater in 1-M than 2-M during the mature phase of the storm after 2–3 h and also exhibit a more unsteady behavior over time (see Figs. 1a,b). The largest vertical velocities overall in both 1-M and 2-M occur during the initial phase. Larger mean updraft velocity also occurs in the convective region in 1-M (see Fig. 12). The convective updrafts in 1-M are stronger than 2-M between 2.5 and 8 km; the downdrafts are slightly stronger between 3 and 6 km.
Differences in the evaporation of rain in the convective region appear to play a key role in differences in the intensity of the convective updrafts between 1-M and 2-M. This is demonstrated with an additional sensitivity test using the two-moment scheme, where the predicted N0r in the convective region (defined as the region within 50 km of the leading edge of precipitation) has a maximum allowable value of 107 m−4 for the calculation of rain evaporation only. If the predicted value exceeds 107 m−4 before this limit is imposed, then λr is recalculated using N0r = 107 m−4 and the given rain mixing ratio for the calculation of rain evaporation. No other microphysical processes are modified. This test allows us to isolate unambiguously the impact of large rain evaporation rates due to the large N0r in the convective region in 2-M. In this simulation, rain evaporation rates are reduced in the convective region due to lower N0r as expected, leading to reduced latent cooling and increased mean convective updraft intensity at midlevels relative to 2-M (Fig. 16). Similarly, Tao et al. (1995) found that neglecting latent cooling associated with the melting of precipitating ice increased updraft strength. Tao et al. suggested that this occurred because of the increase in static stability near and just above the melting layer. Changes in the interaction between the cold pool and ambient wind shear with the reduced latent cooling may also impact the intensity and character of the convective updrafts (Rotunno et al. 1988; Weisman and Rotunno 2004).
The stronger convective updrafts in the sensitivity simulation with reduced evaporation are associated with a reduced rearward flux of positively buoyant air into the stratiform region compared to 2-M at midlevels (Fig. 17a). This reduced flux is associated with a reduction in the strength of the mesoscale updraft in the stratiform region (see Fig. 16). In contrast, there is an increased rearward flux of condensate (Fig. 17b) associated with larger amounts of ice condensate in the convective region (not shown). These results are similar to differences between the baseline 1-M and 2-M simulations described previously.
In the previous section we suggested that many of the differences between 1-M and 2-M could be explained by differences in the rain intercept parameter N0r, which had a large impact on the rain evaporation rates. In 2-M, values of N0r are larger in the convective region and smaller in the stratiform region than the constant N0r = 107 m−4 specified in 1-M (which is the same or similar to values used in many other one-moment schemes for N0r, see Table 1). Note that while most one-moment schemes use a constant value for N0r, this is not an intrinsic feature of these schemes (see discussion in section 2b).
Larger N0r in the convective region compared to the stratiform region near the surface in 2-M is consistent with observations of raindrop size distributions in mesoscale convective systems (note that there is an even greater difference in the modeled N0r between the stratiform and convective regions aloft). Waldvogel (1974) compared drop size distributions between convective and stratiform regions in midlatitude mesoscale rain systems using surface disdrometer measurements, while drop distributions in tropical squall lines have been investigated by Tokay and Short (1996) and Atlas et al. (1999). These studies showed that N0r decreases rapidly during the transition from convective to widespread stratiform rain (e.g., see Fig. 4 in Tokay and Short 1996); Waldvogel (1974) termed this decrease in N0r the “N0 jump.” Tokay and Short (1996) found relatively high evaporation rates in the convective region compared to the stratiform region, and concluded that this resulted from greater numbers of small- to medium-sized drops and fewer large drops in the convective region compared with the stratiform region, for a given rain rate. Ferrier et al. (1995) found that reduction of N0r by a factor of 5 between the convective and stratiform regions, consistent with differences in radar reflectivity for a given rain rate, were needed to simulate a widespread trailing stratiform region associated with a tropical squall line. Here, the 2-M scheme produces a similar distribution of N0r through the microphysical processes, without any tuning, for the reasons suggested by Ferrier et al. (1995): large drop concentrations in the convective region due to collision–coalescence versus smaller drop concentrations in the stratiform region due to the dominance of rain formation by melting of ice. The magnitudes of the modeled N0r at the surface are similar to disdrometer measurements in a continental squall line showing mean values of 2.1 × 106 and 5.5 × 106 m−4 in the stratiform and convective regions, respectively (for fits of the data to exponential distributions; see Figs. 8 and 10 in Uijlenhoet et al. 2003). A more quantitative assessment of the modeled rain size distribution parameters using the 2-M scheme will be performed using recent measurements collected during squall-line events in Oklahoma.
The key point is that no single value of N0r is able to capture the distribution of N0r predicted in 2-M and the N0 jump seen in observations. Thus, tuning the specified N0r in 1-M to values predicted in 2-M for a particular region of the storm produces a bias in the other regions. Both the larger N0r in the convective region and smaller N0r in the stratiform region, relative to the constant value of 107 m−4 in 1-M, appear to be necessary for producing the storm morphology seen in 2-M. This point is illustrated by additional sensitivity tests. A sensitivity test using the 1-M scheme, but with N0r reduced in all regions from 107 to 2 × 106 m−4 (typical of the N0r in 2-M for the trailing stratiform region; see Fig. 9), reduces the rain evaporation rate in both the stratiform and convective regions relative to 1-M (note that further reduction of N0r to 5 × 105 m−4 does not produce a large difference in the results compared with using N0r = 2 × 106 m−4). This run improves several aspects of the simulation; in particular, it produces much more trailing stratiform precipitation than 1-M, as shown in Fig. 18 by the radar reflectivity at t = 6 h (similar results are apparent at other times during the mature phase of the storm). There is also a sharp transition zone of very little precipitation between the stratiform and convective regions, which may be more realistic than 2-M (Biggerstaff and Houze 1993). However, the leading edge of convective precipitation is poorly defined, and there is a narrow maximum of precipitation about 50 km behind the leading edge of the storm associated with a secondary region of convection. The cold pool is also much weaker than in 2-M, and the propagation speed is about 1.5 m s−1 slower during the mature phase. There is a smaller rearward buoyancy and condensate fluxes at midlevels into the stratiform region, consistent with the sensitivity test with reduced evaporation in the convective region that was described in the previous section. Thus, reflectivities are weak above the melting layer compared to 2-M as seen by comparing Figs. 3 and 18.
Conversely, a sensitivity test using the 1-M scheme, but with N0r increased in all regions from 107 to 2 × 108 m−4 (typical of the predicted N0r in 2-M at midlevels in the convective region; see Fig. 9), results in increased rain evaporation and very little surface rainfall in the stratiform region (Fig. 19), despite an increased rearward flux of positively buoyant air into the stratiform region and a stronger mesoscale updraft, relative to 1-M.
These results suggest the importance of capturing the variability of N0r between the stratiform and convective regions. Two-moment schemes allow a more rigorous treatment of N0r (as well as λr) but come at a computational cost due to the additional prognostic variable for number concentration. We emphasize that it may not be necessary to use a 2-M scheme to produce this variability if N0r can be diagnosed in the 1-M scheme in such a way as to capture the N0 jump.
Recent observations analyzed by Uijlenhoet et al. (2003) have also suggested that the shape parameter of the raindrop size spectra, μr, varies between the stratiform, transition, and convective regions of squall lines, with larger μr in the transition region and smaller values in the stratiform and convective regions. However, Uijlenhoet et al. note that the significance of these differences was unclear, given the large statistical fluctuations within each region. Nonetheless, differences in μr can impact the rain evaporation rates. Both the 1-M and 2-M schemes could be modified to account for a variable μr (rather than assuming fixed μr = 0 for the exponential distributions used here); such an analysis is beyond the scope of this paper but should be investigated in future work.
Differences in the treatment of snow and graupel between 1-M and 2-M appear to play much less of a role in the development of trailing stratiform precipitation compared to differences in the treatment of rain. Sensitivity tests using the 1-M scheme, but with two moments predicted for either graupel or snow, do not produce large differences from 1-M. We note that there may be various ice microphysical parameter changes for which this finding would no longer hold true. We also emphasize that the ability of models to simulate trailing stratiform regions may be quite sensitive to other aspects of the ice microphysics, such as the number and type of ice species included in the scheme (e.g., snow, graupel, hail; Fovell and Ogura 1988; McCumber et al. 1991; Ferrier et al. 1995). In general, the sensitivity test using the 1-M scheme, but with prediction of two moments for rain, is able to capture reasonably the results using the full 2-M scheme.
6. Summary and conclusions
A new two-moment microphysics scheme has been implemented into WRF and used to examine the impact of cloud microphysics on the development and evolution of the trailing stratiform region of an idealized 2D squall line. The new scheme is based on the work of MCK and predicts the number concentrations and mixing ratios of five species: cloud droplets, cloud ice, rain, snow, and graupel. Results using the new two-moment scheme were compared to results using a one-moment version of the same scheme. The one-moment scheme assumes a constant size distribution intercept parameter, N0, for each of the precipitation species (i.e., rain, snow, and graupel), similar to most, but not all, current one-moment schemes. In the two-moment scheme, N0 evolves freely from the predicted number concentration and mixing ratio, allowing for a more flexible treatment of the particle size distributions.
While the overall storm morphology was similar between the 1-M and 2-M simulations, the 2-M scheme produced a much more widespread and prominent region of trailing stratiform precipitation, relative to the 1-M scheme. The key factor responsible for this difference was the reduced rain evaporation rate in the stratiform region in 2-M. Larger mean raindrop size (i.e., smaller slope parameter, λr) in 2-M than 1-M in the stratiform region also directly contributed to the larger radar reflectivity. A secondary factor was the increased rain evaporation rates in the convective region at midlevels, which led to a reduction in the intensity of the convective updrafts, an increased flux of positively buoyant air at midlevels from the convective to the stratiform region, and an increased intensity of the mesoscale updraft. Faster rates of ice growth associated with the stronger mesoscale updraft in 2-M than 1-M were able to compensate for a large reduction in the rearward condensate flux from the convective region, so that the amount of ice condensate in the stratiform region was similar between 2-M and 1-M (and actually larger in 2-M than 1-M during the early part of the mature stage of the squall line).
Differences in the rain evaporation rate between 1-M and 2-M were the result of differences in the rain size distribution parameters (the intercept parameter, N0r, and λr) between the runs. Predicted values of N0r in 2-M generally ranged from 105 to 107 m−4 in the stratiform region, and from 107 to 109 m−4 in the convective region, compared with a constant value of 107 m−4 specified in 1-M. Larger values of N0r in the convective region were associated with significant collision–coalescence; in contrast, rain in the stratiform region was primarily produced by melting of snow. This variability of N0r predicted in 2-M between the stratiform and convective regions is consistent with surface disdrometer measurements in midlatitude and tropical systems indicating a sharp decrease in N0r between these two regions, referred to previously as the N0 jump (Waldvogel 1974; Tokay and Short 1996; Atlas et al. 1999). Decreasing the value of N0r specified in 1-M to the mean value predicted in 2-M for the stratiform region (2 × 106 m−4) significantly increased the amount of stratiform rainfall due to decreased rain evaporation but produced a different storm morphology relative to 2-M. In particular, this sensitivity test produced a poorly defined leading edge of convective preciptiation, a weaker cold pool, slower propagation speed, and weak reflectivities in the stratiform region above the melting level, compared to 2-M. Conversely, increasing N0r in 1-M to the value predicted in 2-M for the convective region (2 × 108 m−4) resulted in almost no stratiform rainfall due to the increased rain evaporation. The key point is that no single value of constant N0r in the 1-M scheme was able to reproduce the results of the 2-M scheme. It may be possible to diagnose N0r in the 1-M scheme in such a way as to capture the N0 jump; however, the development and testing of such an approach was beyond the scope of this paper.
Because rain evaporation appears to play a key role in the development of the trailing stratiform region, it is likely to be sensitive to additional parameters in the 2-M scheme. Specifically, results may be sensitive to our inclusion of implicit raindrop breakup, parameterized by limiting the slope parameter for rain, λr, to 20 cm−1 [consistent with Hodson (1986), Srivastava (1987), and others], as well as our assumption of rain size distribution shape (i.e., exponential versus gamma with nonzero shape parameter). Recent work has demonstrated the important of the size distribution shape parameter on rain evaporation (Seifert 2008). The prediction of number concentration for snow and graupel had much less impacton results than the prediction of number concentration for rain. Thus, a scheme with two moments predicted for rain only was able to reproduce results using the full 2-M scheme. However, we emphasize that the number and type of ice species (i.e., snow, graupel, and hail) can have a significant impact on development of the trailing stratiform rain region (Fovell and Ogura 1988; McCumber et al. 1991; Ferrier et al. 1995). Examination of the sensitivity to these aspects of the model will be left for future work.
Several previous studies have examined the development and evolution of squall lines using 2D models (e.g., Fovell and Ogura 1988; Lafore and Moncrieff 1989; McCumber et al. 1991; Ferrier et al. 1995; Tao et al. 1995). The justification is that variability perpendicular to the line tends to be much greater than variability parallel to the line, at least over larger scales. However, despite the mostly 2D structure, the 3D flow in and around convective cells may be important for developing and sustaining squall lines (Takemi 2007). We are currently testing the scheme for both idealized 3D cases and a real 3D squall-line case based on recent observations in the central United States; this study will also entail a more quantitative assessment of the raindrop size distribution predicted by the 2-M scheme. Results will be reported in a future publication.
This work was supported by the National Oceanic and Atmospheric Administration (NOAA) Grant NAOAR4310107 (W. Grabowski, PI), Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) Grant DE-FG02-03ER63539 (J. Curry, PI), the National Aeronautics and Space Administration (NASA) MAP NNG06GBB1G (J. Curry, PI), and by the National Science Foundation (NSF) Science and Technology Center for Multi-Scale Modeling of Atmospheric Processes (CMMAP), managed by Colorado State University under Cooperative Agreement ATM-0425247. We thank E. Brandes for helpful discussions, and W. Grabowski and G. Bryan for helpful discussions and comments on the manuscript. We also thank R. Fovell, J. Milbrandt, and an anonymous reviewer for comments that helped to greatly improve the manuscript.
Corresponding author address: Hugh Morrison, National Center for Atmospheric Research, 3450 Mitchell Ln., Boulder, CO 80307. Email: firstname.lastname@example.org
* The National Center for Atmospheric Research is sponsored by the National Science Foundation.
Buoyancy, shown here in thermal units (i.e., normalized by base-state potential temperature), is calculated with respect to profiles that are horizontally averaged over the domain and includes contributions from temperature, water vapor, and hydrometeor perturbations.