The Canadian CloudSat/Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) Validation Project (C3VP) provided aircraft, surface, and remotely sensed observations of cloud and precipitation characteristics to support improved simulation of cold-season precipitation within weather forecast models and new developments in satellite and radar precipitation retrievals. On 22 January 2007, the C3VP campaign executed an intensive observation period to sample widespread snowfall that occurred as a midlatitude cyclone tracked along the U.S.–Canadian border. Surface air temperature and precipitation measurements, combined with aircraft measurement of hydrometeor content and size distribution, are used to examine various assumptions and parameterizations included within four bulk water microphysics schemes available within the Weather Research and Forecasting Model (WRF).
In a simulation of the 22 January 2007 event, WRF forecasts reproduced the overall precipitation pattern observed during aircraft sampling, allowing for a comparison between C3VP measurements and microphysics scheme assumptions. Single-moment schemes that provide flexibility in size distribution parameters as functions of temperature can represent much of the vertical variability observed in aircraft data, but variability is reduced in an environment where the simulated temperature profile is nearly isothermal. Double-moment prediction of total number concentration may improve the representation of ice crystal aggregation. Inclusion of both temperature dependence on distribution parameters and variability in mass–diameter or diameter–fall speed relationships suggest a means of improving upon single-moment predictions.
Advances in computing resources have promoted a consistent trend toward increasing spatial resolution and level of detail in operational weather forecasting models and their related research applications. Currently, the Weather Research and Forecasting Model (WRF) is used within the operational forecasting and research communities and includes various physical parameterizations for the land surface, boundary layer, radiation, and cloud microphysical processes, often adding new schemes and updating previous versions with each new release (Skamarock et al. 2008). In recent years, updates to WRF have provided additional multimoment bulk water microphysics schemes that simulate the development of several hydrometeor species, their interactions, growth, and precipitation. As of the version 3.1.1 release, microphysics schemes ranged from single-moment prediction of warm rain processes (Kessler 1969) to inclusion of multiple ice categories (Lin et al. 1983; Hong et al. 2004; Thompson et al. 2008) and other predicted moments that provide additional degrees of freedom used to characterize particle size distributions (Milbrandt and Yau 2005; Morrison et al. 2009; Hong et al. 2010).
In addition to differences in their number of predicted moments, bulk water microphysics schemes differ in the assumed shape of particle size distributions, assignment of relevant distribution parameters, relationships between particle mass and diameter, relationships between particle diameter and terminal fall speed, and the number of simulated microphysical processes contributing to the sources and sinks of each hydrometeor category. Many single- or double-moment bulk water microphysics schemes assume a form of the gamma distribution for precipitating hydrometeor classes:
where Nx(D) represents the number concentration (m−1 m−3) of particles of a given hydrometeor class (x) and diameter (D), Nox is referred to as the intercept parameter, μx the shape parameter, and λx the slope parameter. Studies of surface rainfall by Marshall and Palmer (1948) and snowfall by Gunn and Marshall (1958) are often cited as justification to simplify the gamma distribution to the exponential case by setting μx to zero. Subsequent discussions herein will focus on the particle size distribution, mass–diameter, and diameter–fall speed relationships for snowfall.
Simulated microphysical processes, vertical motions, and advection contribute to changes in the hydrometeor mass for a specific category and location. A mass–diameter relationship is used to distribute the total mass among particles of varying diameter. For snowfall, schemes examined herein assume a power law following Locatelli and Hobbs (1974):
where the coefficient am characterizes the particle effective density and bm refers to the particle fractal dimension (Liu 1995; Lin and Colle 2011). Similarly, a power-law relationship characterizes the terminal fall velocity as , where the coefficient and exponent vary with crystal habit and degree of riming (Locatelli and Hobbs 1974). Schemes representing snow with fixed-density spheres (ρs) define and bm = 3 (Lin et al. 1983; Tao et al. 2003; Hong et al. 2004; Morrison et al. 2009) although observations indicate that both am and bm vary with ambient temperature, crystal habit, and degree of riming (Locatelli and Hobbs 1974; Heymsfield et al. 2007).
Bulk water microphysics schemes are characterized by the number of moments (Mn) predicted for the particle size distribution of each species within a grid volume. The moment of a size distribution is the product of the particle diameter raised to the moment order and the particle number concentration, summed over the full range of sizes. In general terms and with application of the gamma size distribution in (1), the resulting moments are
By combining (1) and (2) and integrating over the full range of sizes, the total mass concentration for a population of hydrometeors (M) can be determined as a product of the mass–diameter coefficient and :
Following this terminology, single-moment microphysics schemes predicting the total mass concentration of a population of hydrometeors are providing , while schemes with a second moment often include the prediction of total number concentration (M0; Morrison et al. 2009). Other parameters are required to calculate the remaining size distribution variables and are assigned as constants or diagnosed from other predicted variables. Given assignments or parameterizations for Nos, μs, am, and bm, the distribution slope parameter can be determined based upon the product of the predicted snow mixing ratio (qs) and dry air density (ρa):
Since weather radars are used to monitor the spatial coverage, duration, and intensity of precipitation, radar reflectivity is often calculated from forecast model output to depict the evolution of simulated storms. Given the variety in treatments of snow and ice within various bulk water schemes, care must be taken to ensure that scheme assumptions are implemented in the calculation of radar reflectivity. Smith (1984) provides a methodology to determine the equivalent radar reflectivity factor (Ze) for populations of ice particles sampled by C- to S-band operational weather radar systems where Rayleigh scattering assumptions are valid. Calculation of Ze is the product of the radar reflectivity factor (Z), defined as the sixth moment of the particle size distribution of equivalent ice spheres, and the ratio of dielectric factors for ice and water. Through equivalence of a particle with mass determined by (2) and that of a pure ice sphere with density of ρi, the radar reflectivity factor for a gamma distribution of ice particles is
Schemes assuming a fixed bulk density for snow result in a reflectivity factor that is the product of . Observations of snow and ice crystals often produce values of bm around 2.0 (Locatelli and Hobbs 1974), while am can vary with temperature and degree of riming (Heymsfield et al. 2007). Therefore, simulated radar reflectivity calculated from schemes that provide a nonspherical treatment of snow or prediction of am and bm will depend upon their predictions of higher-order moments around M4.
All microphysics schemes include a prediction or assignment of various parameters to characterize the properties of snowfall and the particle size distribution. The means of assigning snowfall characteristics is quite diverse. For example, the WRF 6-class, single-moment microphysics scheme (WSM6; Hong and Lim 2006) retains a fixed density for snow but assigns the distribution intercept as a function of temperature following Houze et al. (1979). This allows for the Nos parameter to increase with colder temperatures aloft and provides an increase in the number concentration of small particles with altitude. The single-moment Stony Brook scheme (Lin and Colle 2011) also incorporates the Houze et al. (1979) parameterization of Nos(T) but includes a diagnosed riming factor (Ri; Lin et al. 2011) and temperature dependence that allow for variability in the parameters of the mass–diameter and diameter–fall speed relationships. Combined with Nos(T), the inclusion of the riming factor permits variations in density and fall speeds for particles ranging from dry aggregates to graupel. The Morrison double-moment scheme (Morrison et al. 2009) predicts the total number concentration for the snow category rather than relying upon Nos(T). This approach provides for the evolution of size distribution characteristics with altitude without relying upon additional model predictions for parameterization. Another alternative is to base size distributions and predicted moments purely upon observations. The single-moment scheme of Thompson et al. (2008) incorporates relationships of particle size distribution moments and air temperature determined by Field et al. (2005), based upon over 9000 samples of ice spectra in various synoptic environments. The unique size distribution shape of Field et al. (2005) is incorporated with the calculation of higher-order moments based upon observed relationships between the predicted snow mass and air temperature. Along with differences in size distribution and mass–diameter relationships, each scheme includes varying numbers of simulated microphysical processes, their order of operations, means of performing saturation adjustment, and relations between particle diameter and fall speed.
The performance of a given microphysics scheme is often evaluated by comparing model simulations to observed events. Shi et al. (2010) evaluated the National Aeronautics and Space Administration (NASA) Goddard Cumulus Ensemble (GCE) scheme (Lang et al. 2007; Tao et al. 2003; Tao and Simpson 1993) for lake-effect and synoptic-scale snowfall events that occurred from 20 to 22 January 2007 and were sampled during the Canadian CloudSat/Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) Validation Project (C3VP). Although precipitation rates were comparable to observations, comparisons for C- and W-band radar reflectivity suggested a bias toward higher reflectivities aloft (Shi et al. 2010). Molthan et al. (2010) evaluated the synoptic-scale snowfall event examined by Shi et al. (2010) and determined that the GCE scheme might be improved upon further by incorporating flexibility in the size distribution, mass–diameter relationship, and terminal velocity–diameter relationship for snowfall. Some of these modifications were also suggested as future work by Shi et al. (2010). The aforementioned WSM6, Stony Brook, Thompson, and Morrison schemes were investigated in this study to evaluate techniques similar to those suggested by Molthan et al. (2010). The C3VP synoptic-scale snowfall event from 22 January 2007 is examined once more to identify other scheme assumptions, parameterizations, or techniques that can improve upon the representation of snowfall within a high-resolution forecast model.
2. The Canadian CloudSat/CALIPSO Validation Project and the 22 January 2007 snowfall event
The C3VP was a collaboration between the NASA Global Precipitation Measurement (GPM) and the CloudSat/CALIPSO science communities. The C3VP campaign collected in situ and remotely sensed observations of snowfall required to improve and validate satellite and radar retrievals of precipitation at high latitude. Deployments of C3VP assets were managed at the Canadian Centre for Atmospheric Research Experiments (CARE). Observations of snowfall were obtained from remote sensing, surface instrumentation, and the National Research Council Canada (NRC) Convair-580 instrumented aircraft. Herein, the NRC Convair-580 flight track is divided into two vertical profiles: a descending spiral over the CARE site and the ascending departure along a southeastern heading (Fig. 1), with each profile including atmospheric state parameters, hydrometeor content, and particle characteristics. Particle size distributions were provided by Particle Measuring Systems Inc. (PMS) 2D-C and 2D-P imaging probes, coincident with measurement of ice water content from a counterflow virtual impactor (CVI; Twohy et al. 1997), and liquid water content (LWC) from Nevzorov, PMS King, and Rosemount Ice (RICE) probes. Aircraft data were postprocessed to 5-s increments of flight time to allow for coincident measurement of particle size distribution and other state data (Nesbitt 2009). Remote sensors included the Environment Canada C-band, dual-polarimetric radar at King City, Ontario (Hudak et al. 2008), the CloudSat Cloud Profiling Radar (Stephens et al. 2002), and passive microwave observations from the Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) and the Advanced Microwave Sounding Unit-B (AMSU-B) (Petersen et al. 2007). Surface observations at the CARE site included disdrometer measurements of the particle size distribution, particle terminal velocities measured by a Hydrometeor Velocity and Shape Detector (HVSD; Barthazy et al. 2004), digital photographs of crystal structures and degree of riming (Petersen et al. 2007), and sensible weather measurements at the experimental U.S. Climate Reference Network (USCRN) site near Egbert, Ontario.
Widespread snowfall occurred across the C3VP operational domain on 22 January 2007, associated with a surface warm front aligned northeast of a midlatitude cyclone that moved across Lake Erie and Lake Ontario from 0000 to 1200 UTC (Fig. 2). The temperature gradient and upper-level warm air advection associated with the boundary were relatively weak. Reflectivity observed from King City and the S-band, Weather Surveillance Radar 88-Doppler (WSR-88D) at Buffalo, New York, monitored the evolution of precipitation from 0500 to 0800 UTC along and north of the boundary as the cyclone progressed eastward and depict the event in the context of available surface observations and the aircraft flight track (Fig. 2). Vertical cross sections of dual-polarization radar parameters were acquired directly over the CARE site along the 331° azimuth of the King City radar. Below 2 km, horizontally polarized radar reflectivity increased sharply, corresponding with a downward decrease in the differential reflectivity (Fig. 3). Combined, these parameters suggest a transition from smaller, horizontally oriented crystals into irregularly shaped aggregates lacking a preferential orientation. Aircraft probe imagery and surface observations confirmed the sampling of aggregates in the dual-polarization radar observations. Based upon these observations, model simulations should represent hydrometeor types dominated by snowfall with a transition in size distribution between the cloud top and the surface, emphasizing changes that occur near the 2-km level. Consistent cloud cover and weak temperature advection minimized the diurnal cycle of temperature, with surface temperatures averaging −9°C during the period (Fig. 4a). Moderate precipitation rates of 0.3 to 0.5 mm h−1 began around 0200 UTC with a storm total, liquid equivalent precipitation of 2.8 mm measured at the USCRN site at 0800 UTC (Fig. 4b).
3. Experimental design and forecast model evaluation
Herein, C3VP surface and aircraft observations from the 22 January 2007 snowfall event are used to evaluate three bulk water microphysics schemes available in the Advanced Research WRF (WRF-ARW) as of version 3.1.1 (Skamarock et al. 2008), each selected to evaluate techniques or capabilities that might improve the representation of snowfall. The Stony Brook scheme was incorporated into version 3.1.1 for purposes of this study and is now available to the community as of the WRF-ARW version 3.3 release. Selected characteristics of these schemes are summarized in Tables 1 and 2. Previously, Shi et al. (2010) and Molthan et al. (2010) used the WRF-ARW to recreate the 22 January 2007 event and determined that the resulting simulation compared favorably against radar observations. Forecasts were initialized at 1200 UTC 21 January 2007 by replicating the Molthan et al. (2010) configuration with initial and boundary conditions provided by 6-hourly analyses from the National Centers for Environmental Prediction (NCEP) Global Forecast System at 1° spatial resolution (Table 3). Simulations were performed on a triply nested grid at 9-, 3-, and 1-km horizontal grid spacing (Fig. 5). The innermost 3- and 1-km domains were configured to provide nesting feedback to their parent domains and therefore only the results of the 1-km domain are shown here.
Hourly snapshots of model output were examined to identify periods that compared favorably to the conditions observed during the C3VP intensive observation period. Model output valid at 0600 and 0700 UTC 22 January 2007 represented the moderate snowfall observed by C3VP ground and aircraft instrumentation while also providing a good fit to the aircraft flight time and sampling periods of interest. Simulated radar reflectivity and surface weather conditions at 0600 UTC indicated that each model forecast provided a reasonable position for the surface warm front and precipitation coverage along with placement of temperature gradients and surface wind shifts (Fig. 6). Radar sampling and microphysics scheme assumptions contribute to differences in reflectivity characteristics across the region. For example, the radar composite in Fig. 2 represents combined coverage of the King City radar with 0.8° elevation and coverage of the Buffalo radar with 0.5° elevation. These sampling strategies may return echoes from altitudes above the large aggregate layer indicated over the C3VP site in Fig. 3. Scheme assumptions of particle size distribution shape, predicted snowfall, and mass–diameter relationships also contribute to reflectivity differences between the naturally occurring snowfall and simulated size distributions of aggregates. For a given mass content differences in am or bm in (2) and subsequent characterization of in (6) can lead to reflectivity differences of 5 dBZ or more as the simulated particle size distribution is adjusted toward populations with varying arithmetic mean diameter . Despite differences in radar reflectivity that are attributable to scheme assumptions further evaluated herein, the four simulations represent a broad shield of snowfall associated with the synoptic-scale forcing north of a warm frontal boundary and are comparable to actual weather conditions observed and sampled during the C3VP campaign.
Point measurements available at the CARE site also indicate that the temporal evolution of the simulated event compares favorably to surface observations (Fig. 4). Variability in scheme assumptions, the onset of precipitation, hourly intensity, cessation, and storm total accumulation varied in each forecast run and at each forecast hour. The Morrison scheme outpaced observed precipitation accumulations in the first 4 h but was nearly identical to storm total accumulations ending at 0800 UTC because of an increase in precipitation rate at the CARE site between 0600 and 0700 UTC. Hourly rates from the Thompson and Stony Brook schemes compared favorably to observations through 0600 UTC then deviated from observations as brief, heavier precipitation rates occurred in the area. The WSM6 scheme delayed precipitation onset and produced lighter precipitation amounts, leading to a prediction of storm total precipitation that was less than observed. All forecasts continued precipitation for 2–4 h after snowfall ended at the CARE site. This period of excess precipitation led to continued overestimates from the Morrison scheme while the storm total accumulation from the Stony Brook scheme approached the observed total. The Thompson and WSM6 forecasts produced lesser amounts of precipitation, but all forecasts followed the general trend of observations and produced reasonable hourly rates. Despite variability in hourly precipitation rate, differences in mean precipitation between the four schemes and within the sampling polygon were not statistically significant at the α = 0.05 level. Based upon a general agreement in temperature, precipitation, and spatial coverage of precipitation, model forecast evaluations used conditions averaged among the 0600 and 0700 UTC valid times and within a polygon encompassing the C3VP surface observation site and the bulk of the aircraft flight track (Fig. 6).
a. Temperature and relative humidity profiles
Model initial conditions, temperature advection, radiative processes, and latent heating contribute to the simulated temperature profile. Mean temperature profiles and temperature range were calculated for each forecast. Ranges in temperature were comparable for each scheme and for clarity only the range from the Stony Brook simulation is shown (Fig. 7). Mean temperature profiles were warmer than aircraft data at all levels with a difference on the order of 2°C, but the range in simulated air temperatures was within aircraft measurements. None of the forecasts resolved the inversion in the lowest 2 km despite a relatively high concentration of model vertical levels near the surface. Warm air advection likely contributed to the weak temperature lapse rate in the lowest 3 km but no significant differences in warm air advection were noted among the forecast simulations. Given the consistency of the temperature error among the schemes, it is likely a result of a warm bias in the GFS initial and boundary conditions with lesser impacts from microphysics or radiation scheme assumptions.
Aircraft instrumentation also allows for a comparison of vapor profiles. Each scheme handles the sources and sinks of water vapor with slight differences based upon the assumed order of operations, characteristics of crystals relevant to depositional growth processes, and saturation adjustment procedures. Comparisons between aircraft relative humidities and model output were determined by calculating values from simulated water vapor fields and reporting the maximum relative humidity with respect to water or ice at each model vertical level (Fig. 8). The maximum value is reported rather than the mean and range to avoid the inclusion of isolated, subsaturated areas that may confuse comparisons against aircraft data in an environment where saturation with respect to water was frequently observed. Results from the Thompson and Morrison schemes indicate that saturation (supersaturation) with respect to water (ice) was maintained through 5 km, followed by a sharp reduction in saturation levels near cloud top. Saturation levels in the Morrison scheme were lower than the Thompson scheme in the 5–7-km level, coinciding with the development of cloud ice crystals in the Morrison scheme (Fig. 9). The Stony Brook scheme maintained a profile nearly saturated with respect to water in the lowest 3 km, then began a steady decrease with altitude. The WSM6 forecast underestimated saturation with respect to water by 5%–10% in the lowest 3 km but the underestimation increased with height and resulted in a profile with maximum water and ice saturation values far below observations at higher altitudes. For short-term forecasts, errors in the vapor profiles of the Stony Brook and WSM6 simulations may be confined to regions where precipitation is occurring but longer-term forecast integrations may produce errors that advect downstream. Profiles that are subsaturated with respect to ice and water may reduce precipitation and cloud cover due to sublimation and evaporation, affecting radiation transfer in downstream locations with greater impacts at longer model integrations.
Differences in simulated saturation levels are likely attributable to the saturation adjustment processes within each scheme. The Stony Brook scheme incorporates a linear decrease in water supersaturation from 0% at 500 hPa to 50% above the 300-hPa level and uses a mass-weighted average of saturation values over liquid and ice following Lord et al. (1984) and Tao et al. (1989). The linear decrease and mass-weighted approach provides for a smooth decrease in saturation with altitude and temperature. The WSM6 scheme relies upon diffusion processes for saturation with respect to ice, followed by condensation or evaporation of cloud water to restore the profile to saturation following Yau and Austin (1979) and Rutledge and Hobbs (1983). The Morrison et al. (2009) and Thompson et al. (2008) schemes incorporate a saturation adjustment step for liquid species and rely upon diffusion processes to handle saturation adjustment in the presence of ice. For this particular forecast and simulated environment, the Morrison et al. (2009) and Thompson et al. (2008) techniques adequately represent saturation (supersaturation) with respect to water (ice). The approaches of the Thompson and Morrison schemes provide a better fit to the overall profile of saturation with respect to ice or water and therefore provide some guidance to improve other schemes.
b. Hydrometeor profiles
Mean profiles of simulated hydrometeor content are shown in Fig. 9, averaged from 0600 to 0700 UTC among all 1-km grid points within the sampling polygon with no thresholds for minimum mixing ratios. The categories of cloud ice, snow, and graupel were summed to produce a mean vertical profile and range of total ice water content for comparison against aircraft CVI measurements. Mean hydrometeor profiles from each scheme indicate that the simulation was dominated by snow within the lowest levels in agreement with digital photographs of large, lightly rimed aggregates at the surface (Petersen et al. 2007). Small amounts of cloud water may have been present based upon crystal probe imagery and surface photography of lightly rimed aggregates (Petersen et al. 2007; Molthan and Petersen 2011). Only the Morrison scheme produced a mean profile with a trace of cloud water near 3 km, but the Stony Brook and Thompson schemes also produced cloud water of 0.01 gm−3 or greater, scattered throughout portions of the sampling polygon. Cloud water in the Thompson scheme was confined to the 2-km level, while the Stony Brook scheme produced small amounts of cloud water through 5 km. The WSM6 scheme did not produce cloud water at any level. Small amounts of simulated cloud water are reasonable given that crystal imagery and photography suggested light riming, but perhaps from liquid water contents that were not measured by the available instrumentation.
Predicted total ice water contents were comparable in magnitude and profile shape among the four schemes but each partitioned the total ice water content in different ways based upon assumed characteristics and microphysical processes for cloud ice and snow. The Stony Brook and Morrison forecasts generated a layer of cloud ice from 4 to 7 km but rapidly transitioned to snow or precipitating ice near 5 km through continued deposition and collection processes. The Thompson scheme was dominated by snow throughout the cloud depth. Meanwhile, the WSM6 forecast retained a vertical profile of cloud ice and snow mass that increased from cloud top to the 3-km level. The range of simulated total ice water contents was within aircraft observations at altitudes of 2 km and above and increased with altitude. Mean profiles of total ice water content generally agreed with CVI observations at altitudes above 3 km, but mean simulated values of 0.1 to 0.2 g m−3 were less than observations in the lowest 2 km of sampled aircraft data. Uncertainty in CVI estimates of ice water content were reported as 11% (23%) at 0.2 g m−3 (0.01 g m−3) by Heymsfield et al. (2005). Subject to measurement uncertainty of 11%, CVI measurements still exceed the mean and maximum values of total ice water content predicted in the lowest 2 km of each simulation. Vertical profiles of total ice water content were comparable despite the differences in assumptions among the various schemes. Although schemes differ in their characterization of mass among categories of snow and cloud ice, mean profiles of total ice water content were less sensitive to individual scheme assumptions of physical processes and particle size distributions.
c. Size distribution parameters
Given a predicted hydrometeor content, microphysical processes relate to the distribution of mass among particles of varying size. Number concentrations of each size depend upon the presumed size distribution and relevant parameters. Aircraft 2D-C and 2D-P crystal imaging probe data provided measurements of particle size distributions, which were postprocessed by A. Heymsfield (2012, personal communication) to accommodate their overlap in size range and to mitigate false counts of small particles that result from the shattering of larger crystals and aggregates against the probe housings (Heymsfield et al. 2004, 2008). Particle size distributions were provided in 5-s increments of flight time along with other measurements available from campaign datasets with size bin centers ranging from 30 to 25 mm (Nesbitt 2009).
In this study, aircraft-measured particle size distributions were used to estimate values of Nos and λs for comparison to microphysics scheme assumptions, predictions, or diagnosed parameters. Smith (2003) describes a technique for the estimation of exponential size distribution parameters for raindrops based upon the moments of measured particle size distributions. Although naturally occurring particle size distributions may deviate from the exponential fit, estimates of Nos and λs are required for comparison to scheme assumptions. The Smith (2003) technique was used to estimate Nos and λs from values of M2 and M3 acquired from particle size distributions truncated to bin sizes of 90 μm or larger. Although the aircraft particle size distributions were carefully examined to eliminate incorrect counts of small particles, the selection of the 90-μm threshold is a secondary means of avoiding false counts of small crystals or fragments. The 90-μm size threshold and subsequent calculations of radar reflectivity or ice water content provided a reasonable comparison to aircraft measurements and CloudSat observations (Molthan and Petersen 2011). Following the estimation of Nos and λs, values were retained if the resulting fit to the observed particle size distribution produces a coefficient of determination (R2) greater than or equal to 0.8.
The truncation of a particle size distribution impacts the calculation of size distribution parameter fits and related moment calculations. Observations of Nos typically increase with altitude and colder temperature (Houze et al. 1979; Ryan 2000) in response to greater number concentrations of small particles, therefore estimates of Nos are influenced by the elimination of small particles via increased in the minimum size threshold. Similarly, the elimination of number concentrations of small particles contributes to a decrease in λs and an increase in the arithmetic mean size. Since M0 represents the total number concentration, truncation at a given particle size will reduce overall counts but the impact of the size threshold is reduced for higher-order moments because of the influence of Dn within the calculation. In each case, truncation leads to a reduction in the moment calculated from aircraft data. Therefore, when aircraft data indicate that a given scheme underestimates a value of Mn, use of the full particle size distribution from the aircraft data would increase the magnitude of the model error. In cases where the model appears to overestimate aircraft data, use of the complete particle size distribution may improve the relationship between observations and model simulation.
Aircraft vertical profiles of Nos and λs are shown in Figs. 10a,b, based upon estimates from aircraft particle size distributions truncated to particles with maximum dimension of at least 90 μm and the Smith (2003) technique. Estimated values of Nos and λs decrease between cloud top and cloud base coincident with continued depositional growth in an ice-supersaturated environment and the development of aggregates. Radar observations suggesting the presence of an aggregate zone are corroborated by reductions in Nos and λs, indicating that number concentrations of small particles were replaced by larger aggregates (Fig. 3). The lowest 2 km of the descending profile sampled a narrow, short-lived band of heavier precipitation likely dominated by large aggregates, a feature that was isolated, transient, and not well represented in the 0600 or 0700 UTC forecast hours. The remainder of the descending aircraft profile is similar to the departure ascent, suggesting that the profile is representative of conditions in the broader precipitation shield sampled by the ascent profile. To better characterize the change in particle size, the mass-weighted diameter is shown in Fig. 10c, determined as . As in Molthan and Petersen (2011), aircraft estimates of Dm assume bm = 2.0, as this provided a reasonable match between CVI estimates and integrations of ice water content from individual particle size distributions. Mass-weighted mean particle size increased downward because of the effects of aggregation with greater Dm again noted among large aggregates in the slowest 2 km of the descending profile.
Aircraft estimates of Nos, λs, and Dm were compared against model assumptions within each of the selected schemes. The WSM6, Stony Brook, and Morrison schemes assign an exponential size distribution to characterize particle size distributions of snow or precipitating ice. Although the Thompson scheme incorporates the size distribution fits of Field et al. (2005), equivalent estimates for exponential parameters are determined from , , and (G. Thompson 2012, personal communication). Mean profiles and range of all parameters were based upon calculations consistent with scheme assumptions. The Houze et al. (1979) parameterization of Nos(T) used within the WSM6 and Stony Brook schemes provided an increase in Nos with height as temperatures decreased aloft (Fig. 10a). Variability in aircraft estimates of Nos increased in the 4–6-km altitude range but mean profiles of predicted Nos were slightly less than aircraft estimates. The single parameterization of Houze et al. (1979) will not fit all cases given the variability in Nos(T) relationships noted in Ryan (2000) but it generally followed the observed increase in Nos with height. The observed temperature profile included a shallow temperature inversion near 1 km and a nearly isothermal layer from 2 to 2.5 km. The simulated temperature profile did not provide the same level of detail, instead predicting a vertical column with a weak lapse rate (Fig. 7). With little change in the simulated air temperature over the 1–3-km depth, Nos(T) was stagnant while aircraft estimates of Nos continued to decrease toward the surface. In the Thompson scheme, the vertical profile of M2 is obtained from the predicted snow mass and M3 is diagnosed as a function of M2 and temperature. In the lowest 3 km, changes in snow mass and air temperature were reduced, limiting changes in Nos with height. At 4 km and higher, vertical variability in predicted Nos improved, but predicted values of Nos were consistently less than aircraft estimates. The double-moment Morrison scheme determines λs as a function of the predicted total number concentration (N, or M0) and snow mass content, then determines Nos as the product of N and λs (Morrison et al. 2009). Among the four schemes, the double-moment Morrison scheme provided a better representation of Nos in the lowest 3 km by avoiding a dependence upon the simulated temperature profile.
Within the WSM6 and Stony Brook forecasts, profiles of λs are constrained by the diagnosed value of Nos, predicted snow mass content, and assumed mass–diameter relationship (5). In the WSM6 forecast, a large fraction of the total ice water content was categorized as cloud ice, but the parameters of the exponential size distribution apply only to predicted snow mass. Cloud ice simulated by the WSM6 forecast was a large contributor to the total ice mass, but was composed of smaller particles not included in calculations of λs. Mean profiles of simulated λs generally overestimated aircraft estimates, though aircraft estimates fell within the range of simulated values. When the cloud ice and snow categories were combined and assumed to fit a particle size distribution with Nos(T), the mean profile of λs was reduced. Resulting mean values of λs were a better fit to aircraft data in the lowest 2 km, and although the mean profile is less than the bulk of aircraft observations at 2 km and higher, the range in the diagnosed λs was within aircraft estimates (Fig. 10b). The Stony Brook scheme incorporates a riming factor and temperature dependence in assignment of the mass–diameter relationship, and when combined with the Houze et al. (1979) prediction of Nos(T), contributed to a profile of λs that followed the general trend observed in the ascending profile. Variability in am and bm with temperature had the greatest impact as the scheme predicted negligible amounts of riming confined to the lowest 1.5 km (Fig. 11a). The relatively shallow layer of predicted riming implies that variability in the mass–diameter and diameter–fall speed relationships are driven by temperature at altitudes above 2 km (Figs. 11b,c). Saturation with respect to water and ice were underestimated at higher altitudes and likely reduced further increases in the riming factor. Given the improved fit to Nos in the lowest 2 km, the Morrison scheme provided a good representation of continued aggregation represented by further decreases in λs. At higher altitudes, the Morrison scheme mean profile underestimated aircraft estimates of λs, coincident with a general underestimation of Nos and perhaps other impacts resulting from the fixed density, spherical treatment of snow used within the mass–diameter relationship. The Thompson scheme produced a mean profile and range in λs that are less than aircraft observations but an increase in λs with height that followed the overall vertical trend.
Vertical profiles of mass-weighted mean diameter (Dm) incorporate the characteristics of the mass–diameter relationship and size distribution parameters. The WSM6 and Stony Brook schemes produced mean profiles of Dm that follow the downward increase in aircraft-estimated values, each predicting a range that is within the bulk of aircraft estimates (Fig. 10c). Combination of cloud ice and snow and assumption of exponential parameters in the WSM6 scheme led to a mean profile of Dm that exceeded aircraft observations at altitudes of 2 km and greater. For the Thompson scheme, the mean profile and range in Dm underestimated aircraft observations. The Morrison scheme produced a mean profile of Dm larger than aircraft estimates throughout the bulk of the aircraft vertical profile, despite a good fit to aircraft estimates of Nos and λs. The WSM6 and Morrison schemes assume a spherical mass–diameter relationship with bm = 3, resulting in prediction of Dm based upon higher-order moments than in either the Thompson (bm = 2) or Stony Brook schemes. Although the Morrison scheme provided a good representation of Nos and λs, the spherical mass–diameter relationship resulted in an overestimate of the mass-weighted mean particle size throughout the bulk of the profile. A similar result occurred when the cloud ice and snow mass categories of the WSM6 forecast were combined. The use of a nonspherical relationship in the Stony Brook scheme improved the fit to observations and is better aligned with natural observations of snowfall with bm near 2.0 (Heymsfield et al. 2002; Heymsfield 2003). The Thompson scheme also incorporates nonspherical treatment of snow but underestimated Dm based upon moments that are predicted from snow mass and temperature. Since the snow mass content was comparable to CVI estimates, the relationships between M2 and temperature in Field et al. (2005) may differ from this specific event.
d. Size distribution moments
In addition to size distribution parameters, the moments of particle size distributions are necessary for the calculation of various microphysical processes. Selected moments predicted by each parameterization scheme were compared against those obtained from the number concentrations and particle sizes measured by aircraft probes during the descending spiral and departure ascent profiles (Fig. 12). As with calculations of size distribution parameters, moments from aircraft data are calculated based upon merged particle size distributions among the 2D-C and 2D-P probes but restricted to particles with maximum dimensions of at least 90 μm. Analysis includes mean profiles for all vertical levels and the range of values simulated at altitudes of 4 km and below, selected to focus on levels where variability in simulated cloud top was reduced and transitions from cloud ice to snow were completed in the Stony Brook and Morrison schemes.
Comparisons of M0, or total number concentration, are shown in Fig. 12a. In the Thompson scheme, M0 is diagnosed as a function of M2 via the predicted snow mass and air temperature, while calculations for the Stony Brook and WSM6 schemes integrate the exponential particle size distribution after determination of Nos and λs. The Morrison scheme explicitly predicts M0 as the snow number concentration. Schemes vary in their own methods for characterizing small, high-density cloud ice crystals and larger, lower-density ice crystals or aggregates. The Thompson scheme specifies 200 μm as the minimum diameter for snow, while the WSM6, Stony Brook, and Morrison schemes transition cloud ice to snow through autoconversion processes. The autoconversion process from cloud ice to snow is triggered when cloud ice mass exceeds a specified threshold related to the maximum permitted size of cloud ice crystals, but it does not require a corresponding minimum size for the snow category. Since the calculation of M0 represents the total number concentration, truncation to particles larger than 90 μm reduces M0. Truncation at larger particle sizes would lead to a greater reduction as additional counts of small particles are eliminated. The Morrison scheme prediction of M0 provided a good fit to aircraft estimates of M0 with a mean profile and range that is within the bulk of aircraft estimates below 4 km and also represented the general decrease from cloud top to cloud base associated with continued aggregation. The WSM6 and Stony Brook formulations produced mean profiles and ranges that were within aircraft observations between 3 and 5 km but mean profiles of M0 provided limited vertical variability in the lowest 3 km, coincident with reduced variability in the vertical profile of temperature and values of Nos(T). The M0 diagnosed in the Thompson scheme produced a mean profile that increased slightly from cloud top toward cloud base, though the range in values is similar from 1 to 4 km. The particle size distributions of Field et al. (2005) that were used to develop parameterizations of M0 as a function of M2 and air temperature did not include strong variability in M0 with temperature, although higher-order moments and characteristic particle sizes consistently decreased at colder temperatures. Given the limited variability in air temperature and snow mass in the lowest 3 km of the simulated profiles, the resulting diagnosis of M0 by the Field et al. (2005) equations may have contributed to the lack of vertical variability and overestimation of M0 in this specific event. In general, the explicit prediction of total number concentration in the Morrison scheme was helpful in this case by continuing the downward decrease in M0 associated with aggregation without depending upon the predicted vertical profile of air temperature or snow mass content.
In the Thompson scheme, M2 is dependent upon the predicted snow mass and the specified mass–diameter relationship. Since the Thompson scheme provides a reasonable depiction of the mean profile of snow mass content, the mean vertical profile of M2 follows the overall vertical trend in aircraft estimates. Although the WSM6 scheme prediction of snow mass and diagnosed values of M0 were within aircraft estimates, performance for M2 was degraded unless the cloud ice mass was included within the calculation of λs. The mean vertical profile and range of M2 diagnosed from the Stony Brook scheme demonstrated increased vertical variability over diagnosis of M0 owing to the downward increase in snow mass content and variability in λs. Although the Morrison scheme predicted a downward decrease in M0 coincident with aircraft observations, the mean profile and range in M2 became stagnant in the lowest 2 km and resulted in values slightly less than aircraft estimates.
Aircraft estimates of M4 and M6 are shown in Figs. 12c,d and serve as evaluations of higher-order moments related to simulated values of radar reflectivity. As noted previously, the descending aircraft spiral sampled a region of enhanced aggregation in the lowest 3 km, which corresponds to a sharp increase in M4 and M6 due to greater number concentrations of larger aggregates. Schemes evaluated here did not predict this feature in the specific location and time sampled by aircraft, but it is included as a reference to the range of aircraft-observed values. At altitudes of 3 km and higher, the descending spiral produced observations comparable to the broader region of snowfall sampled in the departure ascent. Although the mean profile and range of M4,6 in the Thompson scheme underestimated aircraft observations in the lowest 3 km, values diagnosed by the scheme provided a good representation of the overall trend in M4,6 with height. Comparisons for the Thompson scheme improved at 3 km and higher. This corresponds to the good representation of the mean snow mass profile and Field et al. (2005) parameterizations based upon supporting observations of decreases in higher-order moments as functions of decreasing M2 and colder air temperatures. The diagnosis of M4 and M6 by the WSM6 scheme produced a mean profile and range that underestimated observations in the lowest 3 km and required the inclusion of cloud ice in λs to approach aircraft estimates. Inclusion of cloud ice above 3 km produced a mean profile that that exceeded the bulk of aircraft estimates. The Morrison scheme provided a mean profile that was a better fit to aircraft estimates of M4,6 than either the Thompson or WSM6 formulations, with mean values of M6 that are well aligned with aircraft observations in the lowest 3 km of the aircraft ascending profile. Despite the particularly good fit between the Morrison scheme prediction of M6 and low-level aircraft estimates, simulated radar reflectivity was greater than observed from C- and S-band radars observing the event (Fig. 2). These differences may be attributable to the aforementioned choice of spherical, fixed-density aggregates comprising the snow category and their impact on the selection of am and required to simulate radar reflectivity. The Stony Brook scheme produced a mean profile and range of diagnosed M4,6 that provided the best match to aircraft estimates in the lowest 3 km. Although the mean profile tends to overestimate aircraft estimates above 3 km, the range in simulated values within aircraft estimates. Although the match of M4 provided a good fit in the lowest levels, simulated radar reflectivity was often less than observed from adjacent radars. Similar to the case of the Morrison scheme, values of am and bm influence simulated reflectivity and contribute to observed differences.
e. Terminal fall speed relationships
Simulated precipitation rates depend upon the mass–diameter relationship for a hydrometeor class, the particle size distribution, and the assumed relationship between particle size and terminal velocity. The schemes examined in this study assume a power-law fit between particle size and terminal velocity following Locatelli and Hobbs (1974), although the Thompson scheme incorporates an exponential decay parameter to allow for a decrease in fall speed with increasing size (Table 2). Terminal velocities of aggregates were measured at the surface with a HVSD, which provides an estimate of fall speeds with errors estimated as 10% or less for ice particles (Barthazy et al. 2004). Resulting particle fall speeds were provided by GyuWon Lee (McGill University) with care taken to ensure that measurements were not adversely affected by local winds. It is assumed herein that the reported fall velocities are the terminal velocities for each crystal. Reports for over 11 000 individual flakes were accumulated from 0200 to 0800 UTC 22 January and compiled into a joint histogram based upon HVSD measurement size bin and 5 cm s−1 increments in terminal velocity. A large amount of scatter is present in the data because of variability in precipitation intensity during the observation period, variability in aggregate characteristics, degree of riming, and resulting fall speed (Fig. 13). To compare observations to parameterization scheme assumptions, the median terminal velocity was calculated for each size bin reporting at least 50 observations of sizes 1 mm and larger. Size and count restrictions were implemented to account for limitations in the HVSD detection of small particles and to obtain a reasonable sample size.
The best fit power-law relationship between the HVSD size bin dimension and corresponding median terminal velocity produced values of and bυ = 0.145, similar to values for “unrimed radiating assemblages of dendrites” reported by Locatelli and Hobbs (1974). These similarities can be attributed to the presence of large aggregates in the lowest levels observed by the King City radar, particle imaging probes, and digital photography of aggregates at the surface (Petersen et al. 2007; Molthan et al. 2010). The WSM6 and Morrison schemes assign equivalent values for aυ and bυ and these parameters produce terminal velocities greater than the median observed value for diameters larger than 1.5 mm (Fig. 13). Although individual V(D) fits are compared here, schemes simulate precipitation via mass- or number-weighted mean fall speeds based upon the full particle size distribution simulated within each grid volume. Therefore, departures between assumptions and observations at large particle sizes are offset by the greater number concentrations of smaller particles within the exponential size distribution, or thresholds are applied to reduce fall speeds that are considered too large. The Morrison scheme restricts mass- or number-weighted fall speeds of populations of snow crystals to 1.2 m s−1 or less to reduce the overestimation of fall speeds by populations of large aggregates. The Thompson scheme overestimates fall speeds for particles larger than 1 mm but reduces some error for particles larger than 3 mm through the inclusion of the exponential decay of fall speed with increasing size. Since the Stony Brook scheme accounts for variable aυ and bυ as a function of riming and air temperature, the diameter–fall speed relationship will vary by location and altitude (Figs. 11d,e), although the impact of riming is insignificant for this specific case. When averaged among all points within the sampling polygon and over the two hour forecast period, the Stony Brook relationship produced a fit (, bυ = 0.194) that is comparable to surface observations but with a consistent, positive bias of approximately 0.1 m s−1. Since the V(D) relationships among schemes examined here tend to overestimate fall speeds for the range of particle sizes considered and observed, these relationships could be adjusted to improve the model representation of snowfall simulated during the event. Reducing fall speeds may retain populations of snow crystals within growth environments and encourage increases in total ice water content that were underestimated in the lowest 2 km.
4. Summary and conclusions
Continued advancements in computing allow for the increased use of single- or multimoment bulk water microphysics schemes within operational forecast models. Each scheme includes a varying number of simulated hydrometeor classes, related microphysical processes, and assumptions regarding the size distribution, mass–diameter relationship, and diameter–fall speed relationship for each category. Although model verification of surface weather elements is provided by observing station networks, the verification of cloud properties is often limited to field campaigns such as C3VP where intensive observation periods provide either remotely sensed or in situ measurements of hydrometeor characteristics for comparison to model simulations.
During the C3VP campaign, the 22 January 2007 synoptic-scale snowfall event was sampled from both aircraft and surface instrumentation. The WRF model was configured for the simulation of the event, building upon previous analysis by Shi et al. (2010) and Molthan et al. (2010) by investigating four additional microphysics schemes representing diverse assumptions and predicted moments. Through comparisons of aircraft data and model vertical profiles, the assumptions of the single-moment Thompson, single-moment Stony Brook, single-moment WSM6, and double-moment Morrison schemes were examined. The Thompson and Morrison schemes provided the best representation of the vertical profile of water vapor and resulting levels of saturation with respect to water and ice and could therefore serve as guidance for improving either the Stony Brook or WSM6 schemes, which underestimated aircraft observations. All four schemes generated a similar profile of total ice water content and a fair representation of conditions above 3 km, but underestimated total ice water content in the lowest 2–3 km where large aggregates were observed. Although the Thompson, Stony Brook, and Morrison schemes categorized the vertical column as primarily consisting of snow, the WSM6 forecast simulated a large quantity of smaller cloud ice crystals. When exponential size distributions were fit to snow alone, the WSM6 forecast often underestimated observed values acquired from aircraft, however, combining snow and cloud ice mass into the exponential distribution suggested an improvement. Therefore, the WSM6 scheme could benefit from a greater emphasis on the development of snow rather than relatively large quantities of cloud ice.
Schemes that incorporated temperature dependence for size distribution parameters or moments provided vertical variability similar to observations from aircraft data, although other assumptions of fixed density, spherical snow, or underestimation of ice water content often contributed to an underestimate of exponential size distribution parameters. The prediction of total number concentration by the Morrison scheme may improve the representation of aggregates present at 2 km and below. Although the Thompson scheme incorporates temperature-dependent functions to predict additional moments based upon ice water content, the lack of vertical variability in the lowest 2 km of the simulated temperature and snow mass profiles may have limited success in this specific case, particularly for diagnosis of M0. Surface measurements of particle fall speeds suggest that the diameter–velocity functions of the Thompson, Morrison, and WSM6 schemes overestimated fall speeds for sizes larger than 1 mm, while the inclusion of temperature variability within the Stony Brook scheme produced fall speeds closest to observations.
Although field campaign measurements needed for microphysics evaluations are infrequent, these opportunities can be used to evaluate scheme performance to identify opportunities for improvement, which may in turn lead to improvements in quantitative precipitation forecasts and simulated cloud properties. The onset of simulated precipitation will depend upon the characterization of hydrometeors within the vertical column, the distribution of water or ice mass among number concentrations of particles of varying size, and their resulting fall speed. Although storm total precipitation rates were similar among the schemes investigated here, an analysis of model assumptions against aircraft data demonstrated some differences in the categorization of hydrometeors within the vertical column, means of assigning their size distribution, and fall speed. Through appropriate tuning of model performance against in situ field campaign measurements and additional verification studies in other events or weather regimes, the aforementioned schemes can be improved upon to better represent observed conditions provide improved short-term predictions of cold-season precipitation. In addition, representative profiles of ice water content, density characteristics, and particle size distributions can assist in the development of precipitation retrievals when combined with accurate simulation of corresponding brightness temperatures or other remotely sensed quantities. Future efforts will seek combinations of strengths and elimination of weaknesses in microphysics simulations to advance the prediction of precipitation associated with heavy snowfall events and related support of future radar and satellite retrievals.
Model simulations were performed on the NASA Discover Cluster. Data from the HVSD instrument were provided by GyuWon Lee of McGill University. Aircraft data were provided through S. Nesbitt (University of Illinois) with particle size distributions and probe data quality controlled by A. Heymsfield (NCAR). Prime funding for aircraft studies during the Canadian CloudSat/CALIPSO Validation Project was provided by the Canadian Space Agency. Support for data acquisition and related investigators was provided by the Global Precipitation Measurement Project (via Dr. M. Schwaller) and the NASA Precipitation Measurement Mission (via Dr. R. Kakar). This study was supported in part by the National Science Foundation under Grant ATM-0908288 (Colle). The lead author would like to thank Greg Thompson of NCAR for several comments and helpful guidance that improved the analysis and discussion. Contributions from a second, anonymous reviewer were also helpful for improving the clarity of discussions regarding aircraft data analysis and scheme comparisons.