In this study the Weather Research Forecast model is used with 1-km horizontal grid spacing to investigate the microphysical properties of Arctic mixed-phase stratocumulus. Intensive measurements taken during the Department of Energy Atmospheric Radiation Measurement Program Mixed-Phase Arctic Cloud Experiment (M-PACE) on the North Slope of Alaska, during 9–12 October 2004, are used to verify the microphysical characteristics of the model’s simulation of mixed-phase clouds (MPCs). A series of one- and two-moment bulk microphysical cloud schemes are tested to identify how the treatment of snow and ice affects the maintenance of cloud liquid water at low temperatures. The baseline two-moment simulation results in realistic liquid water paths and in size distributions of snow reasonably similar to observations. With a one-moment simulation for which the size distribution intercept parameter for snow is fixed at values taken from the two-moment simulation, reasonable snow size distributions are again obtained but the cloud liquid water is reduced because the one-moment scheme couples the number concentration to the mixing ratio. The one-moment scheme with the constant snow intercept parameter set to a value typical of midlatitude frontal clouds results in a substantial underprediction of the liquid water path. In the simulations, the number concentration of small ice crystals is found to be underestimated by an order of magnitude. A sensitivity test with the concentration of ice particles larger than 53 μm increased to the observed value results in underprediction of the liquid water path. If ice (not snow) is the primary driver for the depletion of cloud liquid water, then the results of this study suggest that the feedbacks among ice–snow–cloud liquid water may be misrepresented in the model.
Supercooled liquid water has been observed in Arctic clouds since the 1960s (e.g., Witte 1968; Jayaweera and Ohtake 1974; Curry et al. 1996). However, since the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment, it has become increasingly clear that Arctic stratocumulus clouds can maintain liquid water at temperatures as cold as −34°C (e.g., Hobbs and Rangno 1998; Intrieri et al. 2002; Shupe and Intrieri 2004; Zuidema et al. 2005). These mixed-phase Arctic clouds (MPCs) are observed over 40% of the time, with maximum frequency observed during spring and fall (e.g., Shupe et al. 2006). Unfortunately, most weather and climate models fail to produce Arctic clouds that maintain liquid water at low temperatures (e.g., Curry et al. 2000; Inoue et al. 2006; Prenni et al. 2007; Sandvik et al. 2007; Klein et al. 2008). This inability to maintain liquid water in Arctic clouds has a dramatic impact on the surface energy budget, since liquid water in clouds causes an increase in the downwelling longwave radiation and a decrease in the incoming shortwave radiation relative to glaciated clouds (e.g., Curry et al. 2000).
Recently, a two-moment bulk cloud microphysics scheme (i.e., predicting the number concentrations and mixing ratios of the various cloud and precipitation species) implemented in the fifth-generation Pennsylvania State University–National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5) has been shown to reproduce springtime Arctic stratocumulus observed over pack ice during the SHEBA experiment (Morrison and Pinto 2005) and fall stratocumulus observed during the Department of Energy Atmospheric Radiation Measurement Program’s (ARM) Mixed-Phase Arctic Cloud Experiment (M-PACE) on the North Slope of Alaska (NSA; e.g., Verlinde et al. 2007). These studies demonstrate that the two-moment microphysics scheme allows for the maintenance of liquid water in the cloud, resulting in a more realistic simulation of downwelling longwave and shortwave radiation. Morrison and Pinto (2006) demonstrate that, relative to simpler one-moment schemes, their two-moment cloud microphysics schemes caused enhanced radiative warming at the surface and an enhanced precipitation rate. Most of the improvements using the two-moment schemes were related to its more detailed treatment of the snow size spectra, which resulted in a smaller snow particle number concentration relative to that diagnosed in the one-moment schemes. Other studies have shown that reduced ice particle concentrations allowed for the maintenance of liquid water in Arctic MPCs (Pinto 1998; Harrington et al. 1999; Jiang et al. 2000; Morrison et al. 2003; Prenni et al. 2007), although these studies did not focus specifically on the treatment of larger particles (snow).
We utilize a two-moment bulk microphysics scheme in the Weather Research Forecast (WRF) model version 2.2 (Morrison et al. 2009a) to study the microphysical properties of Arctic stratocumulus clouds. Our studies follow simulations of clouds that formed in winter storms during upslope conditions in the Colorado Front Range by Reisner et al. (1998). The Reisner study found that the presence of water in the clouds was sensitive to the treatment of the number concentration of snow, consistent with Morrison and Pinto (2006). An overestimate of the snow number concentration resulted in conversions of water vapor to snow (depositional growth) and cloud water to snow (riming) that limited the amount of water in the cloud. Observations of particle spectra taken during the Winter Icing and Storms Project indicated that the intercept parameter for snow, N0s, is approximately 2 × 105 to 6 × 105 m−4 rather than the 2 × 107 m−4 used in previous studies such as Dudhia (1989). The Reisner study found that microphysical schemes that used a diagnostic relationship between N0s and snow mixing ratio or a prognostic equation for N0s adequately simulated the observed values of N0s and were therefore better able to reproduce the observed cloud water content.
Estimates of N0s differ from case to case and vary in space and time. For example, Houze et al. (1979) found that estimates of N0s in midlatitude Pacific Northwest stratiform clouds associated with frontal systems are dependent on height (though expressed as temperature) and range between 4 × 106 m−4 and 2 × 108 m−4, while near the United Kingdom Field et al. (2005) find values closer to the values estimated by Dudhia (1989). Pinto (1998) found a N0s value in autumnal clouds over the Beaufort Sea of 2.5 × 106 m−4. Therefore, values of N0s are not known a priori. As demonstrated by Reisner et al. (1998), a prognostic equation for snow number concentration can be used to calculate the mean N0s, but it is not clear if adding this degree of freedom results in more realistic feedbacks between snow mixing ratio and number concentration and between snow and cloud water mixing ratios. In the current study we investigate these relationships in simulations of mixed-phase Arctic clouds.
Intensive ground, remote sensing, aircraft, and radiosonde measurements taken during M-PACE, 9–12 October 2004, are used to verify the microphysical characteristics of the model’s simulation of MPCs. We verify the model with measurements taken by aircraft and at the ARM sites during the M-PACE experiment. We choose Barrow and Atqasuk, Alaska, to verify the model at both a coastal site, where the clouds may be sensitive to air–sea temperature and roughness differences, and an inland site removed from the Arctic Ocean, the source of moisture during this observational period. We identify the extent to which the model is able to simulate the maintenance of liquid water in clouds at low temperatures, the role of ice and snow number concentrations in glaciating Arctic clouds, the spatial distribution of clouds across a coastline, and the in-cloud vertical air motions.
In section 2, we detail the M-PACE experiment and the weather that occurred during this chosen period. Section 3 describes the WRF model and the cloud microphysics scheme implemented in the model. In Section 4 we compare the model to observations at Barrow and Atqasuk and then compare the two-moment scheme to a series of one-moment cloud schemes to identify whether a more detailed treatment of liquid water and ice–snow in clouds allows for a better representation of the Arctic boundary layer and surface radiation. In section 5, we summarize the findings of this study.
a. Synoptic-scale features
Weather during M-PACE was characterized by a strengthening high pressure system north of Alaska that caused air to flow from pack ice over the open Beaufort Sea to the NSA. The Moderate Resolution Imaging Spectroradiometer (MODIS) satellite image, provided courtesy of the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC; Fig. 1) shows roll clouds that extend from the pack ice to the NSA. These clouds are aligned closely to the direction of the geostrophic winds. The wavelength of these clouds is 10–15 km (Fig. 4) and the cloud heights varied from 1 to 1.5 km over Barrow (see Verlinde et al. 2007 for further details of synoptic conditions during this period).
b. Cloud properties
Cloud properties were measured during M-PACE by a suite of instruments flown on the University of North Dakota Citation, as well as by ground-based sites at Barrow and Atqasuk (Verlinde et al. 2007). The instruments and analysis used to create the ground-based observational dataset are documented in Shupe et al. (2008a). Briefly, a high-spectral-resolution lidar is used to identify cloud phase. Retrievals from a vertically pointing 35-GHz, millimeter cloud radar are used to derive both ice water content and ice particle effective radius with estimated uncertainties of up to a factor of 2 and about 50%, respectively (Shupe et al. 2005). In addition, the cloud radar is used to estimate the vertical velocity from the distribution of returned radar power as a function of hydrometeor radial velocity in the radar volume, with an estimated uncertainty of no more than 0.2 m s−1 (Shupe et al. 2008b).
Brightness temperature measurements at 23.8 and 31.4 GHz from a microwave radiometer (MWR) provide estimates of total condensed liquid water path (LWP), with an uncertainty of 25 g m−2 (Westwater et al. 2001), and precipitable water vapor. Profiles of liquid water content (LWC) are computed adiabatically from radiosonde profiles and the cloud boundaries observed by radar and lidar, then scaled by the MWR LWP. To the extent that the actual clouds are subadiabatic, the retrieved profile shape of LWC may be moderately incorrect; however, the vertically integrated liquid is correct to within the LWP uncertainty. Radiation was measured by a complete set of shaded and unshaded longwave and shortwave broadband radiometers. Estimates of LWC, ice water content (IWC), longwave (LW), and shortwave (SW) radiation from this dataset are used to verify the WRF simulations at Barrow and Atqasuk.
The observations taken by aircraft have been analyzed and are documented in McFarquhar et al. (2007). The aircraft measured 513 30-s-averaged size distributions of single layer clouds, 71% of which were in mixed phase. These MPCs were typically dominated by liquid drops with precipitating ice below cloud base. Conclusions of the McFarquhar study that are used to verify our WRF simulations of M-PACE are the following:
The size of water droplets increased from the cloud base to the cloud top. This increase in size with height was attributed to condensational growth and there was an indication that collision–coalescence may have occurred at the cloud top to produce some drizzle.
Most of the ice mass is contained in sizes greater than 1 mm (hereafter ice particle “size” refers to its maximum dimension).
The median mass dimension of the ice increased from the cloud top to the cloud base, indicating that the ice crystals grew as they fell through the cloud and may have grown by aggregation below the cloud base.
3. Model description
a. Model setup
The WRF V2.2 model (Skamarock et al. 2005) is used for this study with three nested grids with horizontal grid spacings of 18, 6, and 1 km (Fig. 1). The boundary layer is well resolved in the vertical by including 20 pressure levels below 800 hPa. The model is forced with lateral and surface boundary conditions using the National Centers for Environmental Prediction (NCEP) Environmental Modeling Center Final Data Global Analysis System forecasts every 6 h (more information is available online at http://wwwt.emc.ncep.noaa.gov/gmb/moorthi/gam.html). The radiation, surface layer, land surface, and planetary boundary layer options used in the model runs are described in Table 1. The model is spun up by integrating from 1200 UTC 8 October to 1200 UTC 9 October. The subsequent 24 h (1200 UTC 9 October to 1100 UTC 10 October) are used in the analysis.
b. Microphysical parameterizations
The microphysical cloud scheme used for the baseline run in this model includes two moments for cloud droplets, rain, ice, snow, and graupel. This means that a prognostic equation for mixing ratio and number concentration is integrated for each of the five hydrometeor classes. Morrison et al. (2009a) provide details of the parameterizations used in this microphysical scheme. The results of the two-moment scheme are compared to simpler one-moment versions of the scheme that do not include prognostic equations for the number concentration of rain, snow, and/or graupel, but are otherwise identical to the two-moment scheme.
The microphysical scheme employed in this study assumes that the size distribution of each hydrometeor class has the form of a complete gamma size distribution:
where D is the hydrometeor maximum dimension and Pc, N0, and λ are the spectral index, intercept, and slope, respectively. For solid phase particles it is assumed that Pc = 0. This means the distribution is a simple exponential function with the largest concentrations at the smallest sizes. Both N0 and λ can be written in terms of the mixing ratio, q, and number concentration, N, of each class:
where m = cDd is the assumed mass–diameter relationship. In our model setup, the particles are assumed to be spherical and d is set equal to 3. The constant c is a function of density and in this study c has two discrete values depending on whether the particle is classified as ice or snow. In most one-moment schemes, N0 is specified and λ is subsequently derived from the predicted q and specified Pc, c, and d by rearranging and combining (2) and (3). Then N can be diagnosed from N0, λ, and Pc following (3). Many existing one-moment schemes specify N0 as a constant for each species (Lin et al. 1983; Rutledge and Hobbs 1983; Dudhia 1989; Grabowski 1998), although a few allow N0 to vary as a function of the local temperature and/or mixing ratio for one or more of the species (e.g., Reisner et al. 1998; Thompson et al. 2008). In our single-moment scheme, unless otherwise stated, N0 is specified using values from Reisner et al. (1998) for rain (1 × 107 m−4) and graupel (4 × 106 m−4) and Grabowski (1998) for snow (1 × 107 m−4; see Table 2). Similar values of N0 for snow of 2 × 107 m−4 were used by Rutledge and Hobbs (1983) and Dudhia (1989), while Lin et al. (1983) used a smaller value of 3 × 106 m−4. We explore sensitivity to the choice of N0 in the one-moment scheme in section 4c.
Bulk microphysical models typically have (at least) two ice hydrometeor classes with densities and fall speeds characteristic of large and small particles in order to simulate the dependence of frozen particle densities and fall speeds on size. From hereon we refer to the small particle mode as ice and the larger mode as snow. In the simulations here graupel does not play a significant role. The bulk density of ice is assumed to be 0.5 g cm−3 and the bulk density of snow is assumed to be 0.1 g cm−3.
c. Experiment design
In Table 2 we list all of the experiments discussed in this paper and values for intercept parameters when specified. The baseline simulation is referred to as two moment (2M). A simulation with no ice initiation (NOICE) is listed in the tables for comparison. We compare the 2M simulation with a series of one-moment simulations:
Snow, rain, and graupel intercept parameters specified with values calculated from 24-h averaged 2M fields, referred to as 1M2M.
Snow intercept parameter specified from Grabowski (1998) and all other intercept parameters predicted, referred to as 1MSN.
Snow intercept parameter calculated from 24-h-averaged 2M fields and all other intercept parameters predicted, referred to as 1M2MSN.
Ice intercept parameter specified to produce ice number concentrations similar to aircraft observations during M-PACE and all other intercept parameters predicted, referred to as 1MICE.
LWP, surface downwelling shortwave radiation, and vertical velocity at 890 hPa in the 1-km domain from the 2M simulation for 1- and 24-h averages are plotted in Fig. 2 to show the structure of microphysical, radiation, and dynamical fields. For the 1-h averages, each of these fields show the spatial variability due to roll clouds forming off the marginal sea ice in the northeast and developing parallel to the geostrophic winds. These roll clouds appear to be in good agreement with the observations (cf. Figs. 1 and 2). However, the transition to open cell convection is not seen in Fig. 2 because we are averaging over 1 h or 1 day. Open cell convection is seen in instantaneous snapshots (Fig. 3). By visually comparing the instantaneous LWP in Fig. 3 with the MODIS image in Fig. 1, it is seen that the modeled rolls appear to be transitioning to open cells too rapidly.
Roll clouds do not form in the 6- or 18-km grids, since the resolution in these grids is too coarse to resolve these features (see below). In addition, the 18-km grid produced stratocumulus at Barrow that was too close to the ground, suggesting that the dynamics associated with the rolls is required to obtain the appropriate boundary layer (BL) depth.
To ensure that our model simulations resolved the mesoscale roll cloud structure we ran two additional simulations with the finest mesh extended to the sea ice edge using 500-m and 1-km horizontal grid spacing. The increased resolution and the extension of the domain to the sea ice edge did not change the structure of the roll clouds near the coast and inland over the NSA.
The roll clouds that form in the model 100 km from the domain boundary to the Alaskan coast have an average wavelength of ∼12 km (Fig. 2a and dashed line in Fig. 4). A MODIS satellite image taken at 0000 UTC 10 October 2004 (with 250 × 250 m2 pixel size averaged to 316 × 316 m2 for this analysis) shows roll clouds with similar 8–12-km wavelengths (Fig. 1 and solid line in Fig. 4), where the average power is estimated using a Morlet wavelet spectrum. The two transects that were used in Fig. 4 are shown in Fig. 1. Transect 1 is the line closest to the Alaskan coast and transect 2 is to the east over the open ocean. The modeled planetary boundary layer height in the region of these transects is approximately 1.1–1.3 km. The aspect ratio of these roll clouds (BL height/wavelength) is approximately 10, which is similar to the aspect ratio of 8 found in the observational paper on roll vortices by Hartmann et al. (1997). An aspect ratio of 10 is large compared to high-resolution modeling studies such as Chlond (1992). However, these studies have domains that only extend ∼150 km from the sea ice, where the aspect ratio of the rolls is expected to be smaller than farther downwind. In our runs with the domain extended to the sea ice edge, the aspect ratio of the roll clouds ∼150 km from the sea ice edge is closer to that found in the Chlond study (an aspect ratio close to 4), despite the inability of our model to resolve turbulent eddies driving the initiation and evolution of the rolls near the ice edge.
Our simulation also captures the transition of the rolls to open convective cells (which is not seen distinctly in Fig. 2 because of time averaging) that is apparent from the MODIS image near and just downwind of Barrow (see Fig. 1). The study by Liu et al. (2006) also successfully simulated boundary layer roll clouds using similar resolution (500-m grid spacing). The turbulent kinetic energy in the 2M simulation along transects 1 and 2 are approximately 1.5 m2 s−2 extending up to 0.5 km in the updraft of the rolls (results not shown), similar to the findings of the Liu et al. study. Since the focus of this paper is on the microphysics, we note the reasonable simulation of the roll structure near the coastline that is important for driving the microphysics, but we leave further analysis of the roll dynamics as the subject of a separate paper.
a. Verification with observations at Barrow and Atqasuk, Alaska
Figure 2b shows that the large LWP in the roll clouds reduces the shortwave radiation that reaches the surface to less than 30 W m−2. The roll clouds only extend to approximately 100 km inland, with a decrease in LWP to less than 120 g m−2 at Atqasuk and a resultant increase in shortwave radiation to values greater than 60 W m−2.
In Table 3 we list values of cloud properties and radiation fields at the surface averaged from 1200 UTC 9 October to 1100 UTC 10 October (ordered from lowest to highest LWP) at Barrow for all the experiments listed in Table 2. Table 4 is the same as Table 3 for fields at Atqasuk. These tables show the sensitivity of the surface longwave and shortwave radiation to the liquid water in the clouds. For LWPs less than approximately 150 g m−2, increased LWP increases the surface longwave radiative flux and decreases the shortwave. For LWP greater than approximately 150 g m−2 both shortwave and longwave are less sensitive to the change in cloud liquid water. For example, the LWP at Barrow increases by 52% when all ice is removed from the clouds (NOICE), but the longwave radiation only increases by less than 1% and the shortwave radiation is essentially unchanged.
Comparisons of hourly-averaged ground-based retrievals of LWC at Barrow over the 24-h period 1200 UTC 9 October–1100 UTC 10 October (Fig. 5a) with the model simulation (Fig. 5b) show that the model underestimates the maximum LWC but closely approximates the height of the cloud top and cloud base, as well as the vertical structure of the LWC, consistent with the adiabatic assumption used in constructing the retrievals. The warm bias in the model is because the NCEP fields used to force the model were approximately 1° warmer than observed. The 24-h mean LWC in the model peaks at 1.3 km with a value of 0.34 g m−3 compared to the retrieved LWC that peaks at 1.15 km with a value of 0.41 g m−3 (Fig. 6). The model underestimates the 24-h mean LWP at Barrow by 32 g m−2 (Table 3) and overestimates it by 76 g m−2 at Atqasuk. Both the modeled and retrieved LWC are within a standard deviation of the measurements taken by aircraft near Barrow (Fig. 6, light shading). A small amount of drizzle (less than 0.03 g m−3) forms during periods when maximum LWC exceeds 0.4 g m−3. The 24-h mean droplet number concentration is a maximum of 33 cm−3 at a height of 1.1 km, which is smaller than mean concentration from aircraft measurements but within one standard deviation (e.g., McFarquhar et al. 2007). However, uncertainties in droplet concentration play much less of a role than uncertainties in representation of ice and snow (Morrison et al. 2008). The number concentration of small ice crystals is substantially underpredicted (by about an order of magnitude) relative to aircraft observations. Possible reasons for and consequences of this underprediction are detailed in section 4c.
Comparisons between the hourly-averaged retrieved and modeled IWC (Figs. 5c,d) show that the model has larger variability on this time scale than the retrievals and is a factor of 2 larger than the retrieved 24-h mean (Fig. 6). The model overestimates the 24-h mean IWP at Barrow by 40 g m−2 (more than factor of 2; Table 3). Large fluctuations occur in the model at Barrow, with IWC varying from approximately 0.02 to 0.07 g m−3. Similar fluctuations occur in the hourly-averaged retrievals with smaller amplitude (ΔIWC = 0.01 g m−3 every 2 h). However, larger-amplitude fluctuations occur on shorter time scales in the retrievals (ΔIWC = 0.2 g m−3 every 10 min; Shupe et al. 2008a), which are not seen in the model.
The variability of cloud condensate is closely coupled with oscillations in the vertical velocity. Figure 7a shows that peak LWC occurs during ascending motion while minima occur during descent. In Fig. 7b we compare the standard deviation of vertical velocity at 1 km (below the peak LWC) in 2M (solid) and the retrievals (dash) for averaging periods from 10 s to 1 h calculated with a running mean average. The standard deviations are calculated from 24-h time series with a uniform time step of 10 s (the retrievals have been interpolated to a uniform 10-s time step). The figure clearly shows that the retrievals have larger variability on time scales shorter than 2 min while the model has variability somewhat larger than observed for time scales between 2 and 35 min. Inability of the model to capture variability over time scales <2 min is expected since it cannot resolve turbulent motion, while the retrievals may indeed represent these scales. However, links between spatial scale and temporal fluctuations of the vertical velocity following Taylor’s hypothesis (Taylor 1938) are complicated by anisotropy. In particular, this analysis of the modeled and observed vertical velocity is likely to exclude most of the variability associated with the mesoscale roll dynamics, because the rolls are aligned parallel to the geostrophic wind and therefore do not advect with the mean flow. As shown in Fig. 2c, the largest fluctuations of vertical velocity in the simulation are indeed associated with the rolls. Nonetheless, the similarity of the modeled and observed roll cloud wavelengths (see Fig. 4) indirectly suggests that the model is able to reproduce at least the gross features of the roll dynamics.
To quantify the relationship between the fluctuations in IWP and LWP, Fig. 8 shows the correlations between LWP and IWP at zero lag for the period 1200 UTC 9 October–1100 UTC 10 October for averaging periods from 10 s to 1 h from retrievals (Fig. 8a) and the 2M simulation (Fig. 8b). The model shows a highly significant correlation between IWP and LWP for all averaging periods while this relationship is only seen in the retrievals for averaging periods less than 14 min. The rapid falloff of retrieved correlations at scales between 2 and 20 min (∼1–10-km spatial scale using Taylor’s hypothesis) is not captured by the model, which likely reflects the fact that the model cannot resolve the energy spectra at scales less than the effective resolution of about 5–7 km. However, we again emphasize that this type of analysis cannot capture variability associated with the mesoscale rolls since they are aligned parallel with the mean flow. Overall, we conclude that the model reasonably reproduces the Arctic MPC mesoscale roll structure and many microphysical features, although there are indications that some of the dynamical forcing of the clouds may be on the wrong time and/or spatial scale in the model, mostly due to the inability to resolve smaller- and turbulence-scale features.
At Atqasuk, the 24-h-averaged LWP is a factor of 2 greater than the retrieved value (Table 4), resulting in an overestimate of surface longwave radiation and an underestimate of surface shortwave radiation. However, this overestimate only occurs during the time period from decimal days 9.5–9.75; after this time both the magnitude and variability of hourly-averaged modeled LWP are quite similar to the observations (not shown). As shown in Fig. 2a, the model produces bands of high LWP that extend inland by 100 km or more, one of which happened to be situated directly over Atqasuk for the initial part of the model simulation. While the observational data at Atqasuk do not indicate the presence of this band, the general distribution of modeled LWP in the vicinity of Atqasuk, specifically just to the south, is in much better agreement with the Atqasuk measurements. This agreement, as well as the similarities between Figs. 1 and 2, indicates that the model has some skill at predicting the cloudiness in the regions where roll clouds are advected over the land surface. However, there is some indication that the roll clouds may extend moderately too far inland.
b. Hydrometeor initiation and growth mechanisms in the control simulation
In Fig. 9 we compare the dominant growth and initiation mechanisms for cloud droplets, ice, and snow at Barrow. Cloud droplets are activated in the model using the grid-scale and subgrid vertical motion (Morrison and Pinto 2005) using a lognormal aerosol size distribution to derive the cloud condensation nuclei spectra following Abdul-Razzak and Ghan (2000; see our Fig. 9a). Dry aerosol size distributions assume a bimodal fit to 10 October Met One Handheld Particle Counter measurements of aerosol size distribution and are constrained by measurements of mean condensation nuclei taken by the National Oceanic and Atmospheric Administration (NOAA) Climate Modeling and Diagnostics Laboratory (CMDL). The bimodal scheme is a function of standard deviation, geometric mean, and total number concentration of a smaller and larger aerosol mode. Further details of the aerosol specifications are provided in Morrison et al. (2008). Activation of cloud droplets primarily occurs at cloud base. The activation at the cloud top occurs in regions of low water content. Maximum values near the cloud base occur during updrafts (Fig. 7) with coincident growth of cloud droplets due to condensation (Fig. 9d). The growth of cloud droplets due to condensation is a maximum at cloud base and above in regions of upward motion. The overlay of LWC contours in Fig. 9d shows that these cloud droplets are advected upward and grow (subgrid mixing and radiative cooling also lead to condensation within the cloud layer). This result is consistent with the aircraft observation that found an increase in liquid effective radius from the cloud base to the cloud top (McFarquhar et al. 2007). The evaporation of droplets causes a net decrease in LWC and a sharpening of the LWC gradients at the cloud top. Also seen in Fig. 9d is a decrease in LWC due to evaporation in the downdrafts. There is some evidence of entrainment and evaporation based on aircraft observations of subadiabatic liquid water content near the cloud top, which cannot be explicitly represented in the model because of the relatively coarse grid spacing but appears to be reasonably simulated by the parameterized mixing scheme (see Fig. 6).
Ice is initiated by condensation freezing, deposition, contact freezing, and immersion freezing in the model. Homogeneous freezing of cloud droplets is negligible for temperatures observed during M-PACE. The concentration of ice nuclei acting in deposition and condensation freezing modes is specified from observations during this case using the continuous flow diffusion chamber (CFDC; Prenni et al. 2007). The dominant ice initiation mechanism in the model is immersion freezing, which converts the generally larger cloud droplets found near the cloud top in regions of supersaturation with respect to ice into ice crystals (Fig. 9b). As is seen in Fig. 9b, ice initiation due to immersion freezing increases from 1 L−1 h−1 at 1200 UTC 9 October to 5 L−1 h−1 at 1100 UTC 10 October, all other ice initiation mechanisms do not exceed 0.3 L−1 h−1 during this period.
Snow is initiated by autoconversion of cloud ice to snow, which is parameterized following Harrington et al. (1995; see our Fig. 9c). Ice and snow are initiated at the cloud top, with snow initiated over a deeper layer, and grow as they fall through the cloud, consistent with aircraft observations (McFarquhar et al. 2007) and ground-based analyses (Shupe et al. 2008a). The initiation of ice and snow increases over the 24-h period as the temperature decreases. The dominant growth mechanism for snow is by riming of droplets from the cloud top to the cloud base. Growth of snow due to vapor deposition occurs at a slightly slower rate than riming within the region of water saturation below the LWC maximum. Both mechanisms for the growth of snow reach maximum rates in updrafts.
c. Sensitivity studies
In Fig. 10 we compare the 2M simulation with two one-moment simulations, 1M and 1M2M, and with retrievals at Barrow over the 24-h period 1200 UTC 9 October–1100 UTC 10 October. The run with the snow intercept parameter set equal to 1 × 107 m−4 and all number concentrations for other species predicted (1MSN) is essentially the same as the 1M run, demonstrating that an increase in N0s alone significantly decreases the LWP due to the increased conversion of vapor to snow and cloud droplets to snow (figure not shown).
The 2M and 1M2M runs closely simulate the observed shortwave and longwave radiation at Barrow, with the most notable difference being a −5 to 10 W m−2 bias in the longwave radiation (Figs. 9a,b). Both 2M and 1M2M produce clouds that maintain liquid water in the presence of ice and snow (Fig. 10c). The 2M LWP is within one standard deviation of the retrievals except for the first 6 h. The LWP in the 1M2M is uniformly smaller than the 2M simulation, with a consistent relative increase in shortwave radiation (Fig. 10a). The 1M2M simulation also has larger fluctuations in IWP relative to the 2M simulation (Fig. 10d). These interesting differences between the two- and one-moment schemes are investigated in more detail later in this section. Figure 10d shows that all of the simulations have IWPs that are more variable than the retrievals at the hourly-averaged time scale, consistent with the high correlation between modeled LWP and IWP shown in Fig. 8. The thin dotted lines in the figure show that, just as was seen in the Reisner et al. (1998) study, setting N0s = 1 × 107 m−4 causes the cloud water to be depleted resulting in large biases in shortwave and longwave radiation. The large increase in IWP in 1M is seen in larger-amplitude spikes and an increase in the time mean.
Cloud water droplets grow by condensation (Fig. 9d). Increasing N0s causes an increase in vapor deposition to snow, making less vapor available for the condensational growth of cloud droplets (results not shown). This finding is similar to the results of Reisner et al. (1998) and Morrison and Pinto (2006). In addition, an increase in N0s increases the riming of cloud droplets onto snow. Setting N0s to a value larger than that calculated in the model results in a large decrease in the LWC due to the increased conversion of vapor to snow and cloud droplets to snow (simulations 1M and 1MSN).
As seen in Fig. 10, the one-moment scheme does a good job of reproducing the observed LWP when N0 for snow, rain, and graupel are calculated from 2M fields (simulation 1M2M) relative to the 1M simulation, but generally underestimates the LWP (with consistent changes in the surface radiation fields). To identify to what extent the estimate for snow number concentration used in the one-moment scheme causes the underprediction of the LWP, we plot in Fig. 11 the hourly-averaged snow number concentration and the snow mixing ratio averaged in the cloud (LWC > 0.01 g m−3) for the 2M and 1M2MSN simulations. As is seen in Table 3, the 1M2MSN simulation does a much better job reproducing the observed MPC than the 1M simulation, but still has 25% too much ice and too little water in the clouds relative to the 2M simulation. As is seen in Fig. 11, the time evolution of the 1M2MSN snow number concentration and mixing ratio is quite different from the 2M simulation. Specifically, fixing N0s does not allow the snow number concentration to increase with time as rapidly as in the 2M simulation. In addition, the 1M2MSN snow number concentration and mixing ratio fields are much more variable than the 2M fields, with large oscillations occurring approximately every 4 h. Fixing N0s to a constant value causes the snow number concentration to be a function of the snow mixing ratio only [Eqs. (2) and (3)]. As was shown in Fig. 8, cloud ice and cloud water mixing ratios are correlated in the 2M simulation. Similarly, the cloud water mixing ratio and snow number concentration are correlated in the 1M2MSN simulation (results not shown). This distorts the feedbacks between snow and cloud droplets by causing large negative cloud droplet mixing ratio tendencies due to diffusional growth and conversion of cloud droplets to snow when the cloud droplet mixing ratio is a maximum. These large negative tendencies are essential for maintaining the large oscillations (the “turn-around” mechanism) seen in Fig. 11 that reduce the amount of cloud water relative to the 2M simulation. While 1M2MSN suggests some success with a priori information on N0s, this comparison clearly shows the inability of this type of formulation with constant N0s to capture the temporal evolution of the clouds.
To identify why the 2M simulation is able to reproduce the observed LWP while the 1M simulation is not, we plot the size distributions of cloud droplets and snow in the 2M simulations in Fig. 12. Based on analysis of the observed spectra of snow particles (for sizes larger than 100–1000 μm) by Morrison et al. (2009b, manuscript submitted to J. Atmos. Sci.), N0s is estimated to range between about 2 × 105 and 4 × 106 m−4 during the 10 and 12 October flights depending on different assumptions used in the exponential size distribution fits to the measurements. Thus, the mean N0s of 7.5 × 105 m−4 appears to be reasonable compared to observations, while the value of 1 × 107 m−4 specified in 1M and commonly used in other one-moment schemes is much larger than the observed N0s.
Despite the fairly good agreement between the observed snow size spectra and the modeled values from 2M [cf. Fig. 12 and McFarquhar et al.’s (2007) Fig. 5], the concentration of smaller ice crystals (represented by cloud ice in the model) is much smaller than observed. However, it is important to point out that the observed concentrations of small crystals are much more uncertain than for larger snow particles since the observations only include particles larger than 53 μm, in addition to possible issues related to phase identification and instrument artifacts such as shattering (see discussion in McFarquhar et al. 2007). The underprediction of small ice particle concentrations is consistent with previous MPACE simulations (Morrison et al. 2008; Luo et al. 2008) and appears to result from the discrepancy between the ice nuclei concentration specified in the model acting via the deposition/condensation freezing mechanisms (based on MPACE CFDC observations) and the observed crystal concentration (Fridlind et al. 2007). Uncertainties in the treatment of contact and immersion drop freezing in the model may also contribute to this underprediction.
To identify the sensitivity of the simulation to ice concentration we ran another simulation with N0i specified at 5 × 107 m−4 to produce a total crystal concentration (cloud ice plus snow) for particles larger than 53 μm close to the observed value, and all other species calculated using the two-moment scheme (1MICE). Even though the total number concentration of ice (modeled snow and ice truncated to d > 53 microns in Fig. 13) is much closer to observations (shown with shading in Fig. 13), this run produced a 78% reduction in the mean LWP at Barrow (Table 3). The increase of ice number concentration in this simulation also results in a smaller increase in snow particle concentration relative to the 2M simulation because of the increased autoconversion of particle number from cloud ice to snow. The depletion of liquid water in the cloud is due to the dependence of diffusional growth and riming on the number concentrations of snow and ice. In addition, because snow and ice number concentrations increase throughout the cloud, conversion of vapor to snow and ice also occurs to a significant extent at cloud bottom in 1MICE, not just at the cloud top as in the 2M simulation, thereby diminishing the water vapor that is the source of growth for droplets. This result is consistent with previous studies of the sensitivity of cloud liquid water in MPCs to ice nuclei concentration (e.g., Pinto 1998; Harrington et al. 1999; Jiang et al. 2000; Morrison et al. 2008).
In this study we used the WRF model run with three nests (18-, 6-, and 1-km horizontal grid spacing) to investigate the microphysical properties of Arctic stratocumulus clouds. Roll clouds observed during M-PACE were realistically simulated in the 1-km domain. The modeled rolls had approximately the correct spacing and extended over the near-coastal land similar to satellite observations. The baseline simulation of hourly-averaged cloud top, cloud base, LWC, longwave and shortwave surface downwelling radiation closely approximated the ground-base observations, while the IWC was found to be 2 times larger and more variable compared to observations. However, the snow and ice number concentrations were approximately an order of magnitude smaller than observed. Correlations between IWP and LWP were significant for all averaging periods from 10 s to 1 h in the model, while the observations only show significant correlations for averaging periods less than 14 min, suggesting that the effects of fast processes are underrepresented in the model and those of slower processes are too prominent. This is consistent with the underprediction of the vertical velocity variance over short time scales and the inability of the model to resolve turbulence.
In the model, cloud droplets were activated primarily at the cloud base and grew primarily by condensation while being advected upward in updrafts. Ice was initiated primarily at the cloud top during lifting motions and grew as the ice fell through the cloud. Liquid water was observed to persist through oscillations in the resolved vertical motion while cloud ice was often nearly depleted during descending motions. This general cloud dynamical–microphysical structure is in good agreement with the conceptual model for these clouds developed from observations (Shupe et al. 2008a), albeit at somewhat different time scales.
We found that the model’s simulation of MPCs is sensitive to the snow intercept parameter, N0s, similar to the studies of Reisner et al. (1998) and Morrison and Pinto (2006). Typical values of N0s taken from midlatitude observations and commonly used in one-moment schemes caused the conversion of cloud droplet to snow to occur too rapidly, depleting the cloud liquid water. Fixing the snow intercept parameter to a value taken from the two-moment integration resulted in a mean LWP that was 20% smaller than the two-moment simulation because of a coupling between the snow number concentration and mixing ratio that resulted in large oscillations in the snow mixing ratio that reduced the mean cloud liquid water. Hence, even if a priori knowledge of N0s were available, a single-moment microphysics scheme using constant N0s would likely produce similar oscillations. It may be possible to capture some of the variability of N0s in a one-moment scheme by diagnostically relating this parameter to one or more model variables (e.g., the snow mixing ratio), although development and testing of such a parameterization is beyond the scope of this study. We also note that tuning of the specified N0s in one-moment schemes to improve results for this particular case may degrade results for other case studies.
The two-moment scheme produced mean snow size spectra that were similar to observations, although N0s was somewhat too small. However, our sensitivity studies demonstrate that the liquid water is dependent upon assumptions made about both the snow and ice intercept parameters (or alternatively the number concentrations). While the snow size distribution can be reasonably constrained by aircraft data, the distribution of smaller ice particles is much less certain from measurements. Nonetheless, a simulation with N0i specified to produce total ice (cloud ice plus snow) number concentration for particles larger than 53 μm similar to observations resulted in a depletion of the liquid water by 78%. Potential ice nucleation mechanisms that could explain the discrepancy between the observed ice nuclei and ice crystal concentrations, assuming no significant errors in the measurements of ice crystal concentration, are discussed by Fridlind et al. (2007). For example, the Fridlind et al. study suggests that the formation of ice nuclei from drop evaporation residuals and drop freezing during evaporation may be two mechanisms that need to be included in models in order to accurately simulate observed ice number concentrations in mixed-phase clouds.
The substantial underprediction of LWP and overprediction of IWP in the simulation with the ice number concentration constrained to be similar to observations indirectly suggests that the growth of ice in the model may be unrealistic, assuming that there are no large, undiagnosed errors in the measurements that would lead to observed crystal concentrations significantly larger than in reality. Possible explanations include uncertainties in the treatment of capacitance and vapor diffusional growth, as well as crystal shape effects and representation of the shape of the particle size distribution. If ice (not snow) is the driver for depletion of cloud liquid water, then the results of this study suggest that the feedbacks between ice–snow–cloud liquid water may be misrepresented in the model.
This research was supported by the Office of Science (BER), U.S. Department of Energy, Grant DE-FG02-05ER63965 and the National Science Foundation Grants ARC0612428 and ARC0714053. We appreciate the insightful comments of three anonymous reviewers.
Corresponding author address: Amy Solomon, NOAA/Earth System Research Laboratory, Physical Sciences Division, R/PSD1, 325 Broadway, Boulder, CO 80305-3328. Email: email@example.com
* The National Center for Atmospheric Research is sponsored by the National Science Foundation.