The OWLeS IOP2b lake-effect case is simulated using the Weather Research and Forecasting (WRF) Model with a horizontal grid spacing of 148 m (WRF-LES mode). The dynamics and microphysics of the simulated high-resolution snowband and a coarser-resolution band from the parent nest (1.33-km horizontal grid spacing) are compared to radar and aircraft observations. The Ice Spheroids Habit Model with Aspect-ratio Evolution (ISHMAEL) microphysics is used, which predicts the evolution of ice particle properties including shape, maximum diameter, density, and fall speed. The microphysical changes within the band that occur when going from 1.33-km to 148-m grid spacing are explored. Improved representation of the dynamics at higher resolution leads to a better representation of the microphysics of the snowband compared to radar and aircraft observations. Stronger updrafts in the high-resolution grid produce higher ice number concentrations and produce ice particles that are more heavily rimed and thus more spherical, smaller (in terms of mean maximum diameter), and faster falling. These changes to the ice particle properties in the high-resolution grid limit the production of aggregates and improve reflectivity compared to observations. Graupel, observed in the band at the surface, is simulated in the strongest convective updrafts, but only at the higher resolution. Ultimately, the duration of heavy precipitation just onshore from the collapse of convection is better predicted in the high-resolution domain compared to surface and radar observations.
Lake-aligned, lake-effect snowbands, such as those observed over Lake Erie and Lake Ontario, often produce heavy snowfall downwind of the lake (Niziol et al. 1995; Veals and Steenburgh 2015). Snow accumulations may exceed 1.5 m in 48 h, with mere centimeters accumulating just a few kilometers away (Buffalo NY Weather Forecast Office 2019b). These significant snow accumulations cause power outages, collapsed roofs, impassible roads and stranded motorists and ultimately loss of life (Schmidlin 1993; Buffalo NY Weather Forecast Office 2019b). These bands are generally narrow (~20 km wide), and therefore, while band-scale circulations can be resolved by models with O(1) km horizontal grid spacings (Steenburgh and Campbell 2017; Campbell and Steenburgh 2017; Bergmaier et al. 2017), finer-scale features, in particular embedded convection, cannot. Lake-effect snowband dynamics and the dynamical impacts on microphysics and precipitation are therefore not entirely represented in models that employ O(1) km grid spacings. Knowledge gained from finescale modeling validated against observational studies of these localized events can be used to improve forecasts for society’s benefit.
Both observations and numerical models have been employed to improve lake-effect snow forecasts, with a focus on heavy-precipitation events. Lasting, heavy snowfall tends to occur when winds blow parallel to the long axis of an elliptically shaped lake (e.g., Lakes Erie and Ontario) (Holroyd 1971; Kristovich and Steve 1995; Niziol et al. 1995), which allows the lower atmosphere to humidify and destabilize and favors strong convergence along the length axis of the lake (McVehil and Peace 1965; Peace and Sykes 1966; Passarelli and Braham 1981; Niziol et al. 1995). These snowbands are referred to as long-lake-axis-parallel (LLAP; Steiger et al. 2013) bands. Owing to the specific synoptic conditions that LLAP bands require and the band-scale circulations needed for maintenance, prior work has focused on understanding how synoptic conditions impact LLAP snowband formation (Holroyd 1971; Lavoie 1972; Hjelmfelt 1990) and the mechanisms by which snowbands form and evolve on meso-α–β scales. Once formed, band intensity and maintenance are tied to the nature of the cross-band circulations.
Observations have been used to analyze the formation and maintenance of LLAP bands including the importance of low-level convergence from both the land breeze (from opposite shores) and friction (McVehil and Peace 1965; Peace and Sykes 1966; Passarelli and Braham 1981). Additionally, observations have highlighted the increase in precipitation when snowbands interact with orography (Peace and Sykes 1966). More recently, coordinated observations during the Lake Ontario Winter Storms (LOWS) project (Reinking et al. 1993) and the Ontario Winter Lake-effect Systems (OWLeS) field campaign (Kristovich et al. 2017) have provided the datasets necessary to use in conjunction with models to improve the prediction of severe lake-effect snow. For example, Veals and Steenburgh (2015) used operational radar observations to conclude that LLAP bands are the second most common mode of precipitation over the Tug Hill located downwind of Lake Ontario, behind “broad coverage.” Minder et al. (2015) used both Micro Rain Radars and an X-band profiling radar to determine that precipitation transitioned from convective to stratiform across the shoreline and farther inland during the OWLeS IOP2b case (a 24-h LLAP-band event during which 103 cm of snow fell on Tug Hill). Their conclusion, that the ~400 m high Tug Hill terrain did not invigorate the lake-effect convection onshore, was corroborated by Welsh et al. (2016) who used aircraft observations along with other measurements to study characteristics of the ~3 km deep IOP2b LLAP band over the lake and onshore. Welsh et al. (2016) found updrafts of up to 10 m s−1 over the lake accompanied by heavily rimed snow. They also noted that convection collapsed onshore and in-cloud turbulence decreased, but snow particle sizes increased from stratiform ascent over the terrain and over a cold dome associated with the polar air mass that came around the south side of Lake Ontario and produced a land-breeze front stretching from the southeast corner of the lake to the northeast, over Tug Hill. Campbell et al. (2016) looked at the impact of terrain during the IOP2b case and found that during times when the banding was weak and more broad in coverage, the orographic enhancement was strongest.
These observational studies motivated high-resolution, O(1) km, modeling studies of the OWLeS IOP2b case to explore the impacts of Tug Hill on LLAP band snow accumulation (Campbell and Steenburgh 2017) and impacts of the mesoscale forcing mechanisms, in particular the shoreline geometry (Steenburgh and Campbell 2017) and the secondary solenoidal circulation in the LLAP band (Bergmaier et al. 2017). Steenburgh and Campbell (2017) used simulations of the IOP2b case to determine that the shoreline geometry along the southeast part of Lake Ontario produced a critical land-breeze front, consistent with Welsh et al. (2016). Campbell and Steenburgh (2017) used the same simulations to discern that this land-breeze front interacted with Tug Hill to produce increased precipitation mainly through deposition and accretional growth of ice particles lifted over the density current and the hill. Bergmaier et al. (2017) found that these O(1) km simulations capture the cross-band circulation very well, as confirmed by airborne vertical-plane dual-Doppler radar data collected along flight legs across the IOP2b LLAP band. They found this circulation to be driven by solenoidal forcing along a shallow land-breeze boundary and enhanced by latent heat released within the band’s updraft over the lake. Over land, the LLAP band circulation weakened but the land-breeze front had its own shallower asymmetric circulation. Additional modeling studies of other lake-effect snow cases have looked at the impacts of lake ice coverage on the location and intensity of precipitation (Wright et al. 2013), the impacts of different boundary layer and surface layer schemes on precipitation (Conrick et al. 2015; Fujisaki-Manome et al. 2017; Minder et al. 2020), the sensitivity of precipitation characteristics to various microphysics schemes (Reeves and Dawson 2013), and the impact of regional data assimilation on lake-effect snow forecasts (Saslo and Greybush 2017).
The above studies have provided insight into lake-effect storms and the controls on precipitation location and intensity, but all of these modeling studies used O(1) km horizontal grid spacings. Thus, the interaction between the snowband dynamics and microphysics and how precipitation characteristics change at subkilometer grid spacings has not been explored. Airborne cloud radar data shown in Welsh et al. (2016) and Bergmaier et al. (2017) at a resolution of O(10) m reveal ubiquitous convective updrafts over the lake (both within the LLAP band and laterally, under the LLAP band anvil) at scales that cannot be captured by a 1-km, convection-permitting simulation. The issue at hand is relevant to our understanding of gray-zone dynamics, when updrafts are not fully resolved, and how this impacts the microphysics.
Additionally, there are questions as to how improvements to the dynamical representation of a simulated lake-effect snowband (though increased model resolution) impact its microphysical evolution and ultimately the distribution and intensity of lake-effect snowfall, and how the LLAP band forms and evolves in a high-resolution simulation without a planetary boundary layer (PBL) parameterization. Finally, moving to higher resolution modeling will enable the exploration of the lesser known aspects of lake-effect systems, including the interaction between shallow helical-roll convection, which first forms near the upwind shore, and the deeper lake-scale LLAP circulation, which forms farther downwind.
In this study, the OWLeS IOP2b case is simulated using a nested model configuration with a 148-m horizontal grid spacing inner domain. Model output is compared to radar and aircraft observations. Differences in the microphysics and dynamics of the band between the 148-m domain and its parent domain (1.33-km horizontal grid spacing) are explored. In section 2 we describe the IOP2b event, the observations, and the model. Comparisons of model output to the observations are discussed in section 3, and the impact of finescale model dynamics on the microphysical evolution of the band are discussed in section 4. A summary and conclusions are discussed in section 5.
2. The OWLeS IOP2b lake-effect case, observations, and simulation setup
a. Synoptic overview
At 0000 UTC 11 December 2013, an intense LLAP band formed and heavy snow fell for the next 24 h (Steenburgh and Campbell 2017). The snowband was associated with a persistent upper-level trough bringing Arctic air to the region. By 1800 UTC 11 December 2013, North American Regional Reanalysis (NARR) data showed that 850-mb temperatures were −16°C over Lake Ontario and the 850-mb geostrophic wind was parallel to the long axis of the lake (Fig. 1). A total of 112 cm of snow fell in Redfield, New York, in 24 h (Buffalo NY Weather Forecast Office 2019a). A more complete synoptic description of the event can be found in Campbell et al. (2016).
Surface observations used in this study include precipitation measurements from the snow study stations at Sandy Creek (dataset produced by Steenburgh et al. 2014b) and North Redfield (dataset produced by Steenburgh et al. 2014a, see Fig. 2b). Both automated gauge and manual snow board measurements were made at both locations. Automated observations were made using a heated Noah ETI weighing precipitation gauge with a single Alter-style windshield and manual observations were made using a Judd Communications snow depth sensor. The snow depth sensor board was wiped clean every 6 h (0000, 0600, 1200, and 1800 UTC). Additionally, a Parsivel disdrometer was located at Sandy Creek and measured the number of particles and their fall speeds in 32 size bins (dataset produced by Phillips and Knupp 2014).
Radar observations used in this study were obtained from the KTYX NEXRAD S-band (10 cm) radar located in Montague, New York (Fig. 2b). The KTYX radar variables used are the Level II base reflectivity (NOAA/NWS/ROC 1991). The radar quantitative precipitation estimation (QPE) is derived from the Level II base reflectivity following (Campbell et al. 2016). All plots using KTYX observations are made using Py-ART (Helmus and Collis 2016).
Airborne in situ and remote sensing observations utilized in this study were obtained from instrumentation aboard the University of Wyoming King Air (UWKA) research aircraft. Rodi (2011) and Wang et al. (2012) discuss the measurement capabilities of the UWKA in great depth. Flight level in situ measurements include temperature, air pressure, humidity, the 3D wind components, cloud and precipitation particle size distributions, and precipitation particle 2D imaging. In terms of remote sensing, the multiantenna W-band (3-mm wavelength) Wyoming Cloud Radar (WCR) was mounted aboard the UWKA during OWLeS and provided vertical cross sections of radar reflectivity and Doppler velocity along flight tracks. W-band radars provide very high spatial resolution, but their signal attenuates in the presence of much cloud liquid water, and scattering by hydrometeors larger than ~1 mm falls in the Mie regime (Kollias et al. 2007). This study utilizes WCR measurements from only the near-zenith and near-nadir beams (hereafter the “up” and “down” beams, respectively), which were obtained quasi simultaneously and sampled at 20 Hz along the flight track and every 15 m along the beams. The maximum unambiguous Doppler velocity is ±15.8 m s−1 and the minimum detectable signal for both beams at a range of 1 km is about −33 dBZ (Wang et al. 2012).
The aircraft motion was removed from the WCR Doppler velocity measurements when UWKA attitude angles caused the beams to deviate from the vertical. As in Bergmaier et al. (2017), the wind profile from a nearby sounding taken at 2018 UTC from Oswego, New York, was also used to further correct for horizontal wind contamination arising from off-vertical antenna pointing angles, due to these attitude deviations. Following these corrections, the resulting velocity measurements from the two beams provide the profile of hydrometeor vertical velocity above and below the aircraft. For a more thorough description of how WCR velocities were processed along flight tracks across the LLAP band during OWLeS, see Welsh et al. (2016) and Bergmaier et al. (2017).
Particle size distributions were measured aboard the UWKA by two optical array probes during OWLeS, each sorting particles into 101 bins of equal width: a Cloud Imaging Probe (CIP, sizing 0.01–2.51 mm, 25-μm bin width) and a 2D-Precipitation probe (2D-P, sizing 0.1–20.1 mm, 0.2-mm bin width). CIP’s first four size bins are ignored because they lack reliability, thus the minimum size used here is 0.1 mm. The size of snow particles is complex; particle probe sizes mentioned in this study refer to the maximum 2D dimension. In situ liquid water content (LWC) was estimated by several probes on the UWKA during OWLeS, and their measurements generally agreed well. Here we use LWC from the Cloud Droplet Probe (CDP), which resolved the size distribution of droplets in the 2–50 μm size range.
c. WRF simulation setup
This study uses the Advanced Research Weather Research and Forecasting (WRF) Model, version 3.8.1 (Skamarock et al. 2008), run on the Cheyenne high-performance computer (Computational and Information Systems Laboratory 2007). The WRF Model dynamics are fully compressible and nonhydrostatic, and the prognostic equations are time-integrated using a third-order Runge–Kutta method. Four domains are used in a nested configuration (Fig. 2a) with horizontal grid spacings of 12 km for domain 1 (d01), 4 km for domain 2 (d02), 1.33 km for domain 3 (d03), and 148 m for domain 4 (d04). A total of 79 vertical levels are used on a stretched vertical grid with 29 vertical levels in the first 3 km AGL. The advective time steps used are 54, 18, 6, and 1 s for each of the domains. Model output is written every 10 min for d03 and every 5 min for d04. The domain top is at 50 hPa and is rigid with a Rayleigh damping layer applied to the uppermost 5 km. The North American Mesoscale Forecast System (NAM) analyses (NOAA/NCEP/NWS/U.S. Department of Commerce 2015) are used for initial and lateral boundary conditions. All four domains are initialized at 1200 UTC 10 December 2013, and the analysis period, following a 12-h spinup, is from 0000 to 2200 UTC 11 December 2013.
The model setup for domains 1–3 is similar to the setup used by Campbell and Steenburgh (2017). Terrain in domains 1–3 is obtained using the U.S. Geological Survey (USGS) 30-s data (USGS 1996). Additionally, these domains use the YSU PBL scheme (Hong et al. 2006) and first-order turbulence closure. The Kain–Fritsch cumulus parameterization (Kain 2004) is used on domain 1 and calculates moisture tendencies for cloud water and rain but not ice. Domain 4 is setup to run in large-eddy simulation (LES) mode in WRF. In LES mode, a 1.5-order turbulence closure scheme is used, the PBL scheme is turned off and diffusion is used for vertical mixing. Domain 4 (d04) will be referred to as the WRF-LES domain. This domain uses the Shuttle Radar Topography Mission (SRTM; Farr et al. 2007) 3-arc-s (90 m) terrain dataset (Fig. 2b). Random potential-temperature perturbations with maximum magnitude of 0.01 K are added to the 36 westernmost domain 4 tiles and up to 700 mb to reduce the turbulence spinup time on the upwind boundary (Mirocha et al. 2014). Analysis of energy spectra (not shown) reveals that turbulence is spun up by the time PBL air is advected halfway across the lake. Note that two-way feedback between domains is turned off, though a sensitivity study (not shown) with feedback on shows little difference in the evolution of the WRF-LES domain snowband. Turning off two-way feedback allows for cleaner testing of sensitivity to the grid spacings in domains 3 and 4.
All four domains use the unified Noah land surface model (Chen et al. 1996), the revised MM5 Monin–Obukhov surface layer parameterization (Jiménez et al. 2012), the Rapid Radiative Transfer Model for GCMs (RRTMG; Iacono et al. 2008), and Jensen ISHMAEL microphysics.1 The Jensen Ice-Spheroids Habit Model with Aspect-ratio Evolution (ISHMAEL) bulk microphysics scheme (Jensen and Harrington 2015; Jensen et al. 2017, 2018a,b) predicts the evolution of ice particle aspect ratio [Eq. (A2)], or shape, for two ice species, planar-nucleated (ice-one) and columnar-nucleated (ice-two) particles. Additionally, a third ice species, aggregates, is predicted.
In ISHMAEL microphysics, ice particle properties including mass, number, shape, density, maximum dimension, and fall speed are predicted and consistently updated by processes such as vapor growth and riming. These particle properties are used to diagnose particle types, for example, dendrites have low densities and aspect ratios, whereas rimed particles and graupel are more isometric (aspect ratios closer to 1). Thus, improvements to the dynamical representation of the simulated lake-effect snowband in the WRF-LES domain can be directly linked to changes in microphysical processes (e.g., nucleation, vapor growth, riming, aggregation) and how these changes impact ice particle properties. Both shortwave and longwave radiation are coupled to the microphysics. A more detailed summary of the microphysics scheme is in the appendix.
Lake-surface temperatures for all domains are updated every 6 h from the Great Lakes Environmental Research Laboratory (GLERL) Great Lakes Coastal Forecasting System (GLCFS; NOAA/GLERL 2013) following Campbell and Steenburgh (2017). Additionally, ice cover is set manually based on information from the GLERL ice-cover analysis for domains 1–3 following Campbell and Steenburgh (2017). Lake Ontario was nearly ice free during the period of interest, though a few grid boxes in domain 3 are specified as ice-covered in Prince Edward Bay, Chaumont Bay, and Henderson Harbor. For simplicity, no lake ice is specified in the WRF-LES domain. We do not expect this to impact our results considering that the ice coverage in those bays is sparse, Prince Edward Bay is near the boundary of the WRF-LES domain, and the band is generally south of those bays during the analysis time.
3. Observation and WRF Model output comparison
WRF Model output from both the 1.33-km domain (d03) and the 148-m domain (d04, WRF-LES domain) is compared to surface, radar and aircraft observations. Radar observations are used to evaluate the band-scale to convective-scale elements in the snowband. Aircraft observations are used to explore the finer-scale dynamics and microphysics within the band. Model-observation comparisons are used to determine the dynamical and microphysical impacts that occur when smaller scales are resolved in this lake-effect snowband simulation. Surface observations are used to validate the radar-derived quantitative precipitation estimation (QPE) and are compared to model output to relate the dynamical and microphysical impacts that occur when relatively smaller scales are resolved to differences in ice particle properties and precipitation at the surface.
a. Band-scale evaluation using KTYX radar observations
The radar-derived liquid equivalent quantitative precipitation estimation (QPE) is diagnosed from a power-law relationship with the base (Level II) reflectivity. The same equation that was used by Campbell et al. (2016), from Vasiloff (2001), is used here, where Ze = 75S2 and where Ze is the reflectivity factor in mm6 m−3 and S is the liquid equivalent precipitation rate in mm h−1. This equation is verified by both manual and gauge surface observations at two locations. The skill of the radar QPE is unknown at locations away from the surface observations. Radar beam overshoot (Cocks et al. 2016), variation in frozen precipitation type (i.e., snow/graupel) and sublimation/growth of ice below the radar beam can all lead to errors in radar-derived QPE. Nevertheless, radar-derived QPE provides a way to study the spatial variability of precipitation.
The manual snow depth sensor observations totals (Table 1) include some KTYX reflectivity-based disaggregation needed to distribute 6-hourly observations into the 22-h storm period (Campbell et al. 2016). These manual observations are considered to be the best estimates of precipitation at Sandy Creek and North Redfield. Level II radar QPE is also verified with gauge data, but undercatch (Rasmussen et al. 2012) was an issue at both North Redfield and Sandy Creek, though it was less of an issue at Sandy Creek when verified against manual board observations (Campbell et al. 2016). The 22-h radar QPE matches manual observations at Sandy Creek, but it underestimates precipitation at North Redfield by 17% (Table 1, see also Campbell et al. 2016).
The WRF-LES domain produces slightly better precipitation amounts at Sandy Creek (27 mm) compared to d03 (24 mm) when verified against manual (32 mm), gauge (28 mm), and radar (32 mm) observations. The WRF-LES domain produces 54 mm at North Redfield, in better agreement with the radar QPE (48 mm) and gauge observations (43 mm) compared to d03 (57 mm) but slightly worse when compared to the manual observations (58 mm).
Both d03 (Fig. 3a) and the WRF-LES domain (Fig. 3b) capture the total accumulated precipitation over the entire 22 h period (also, see Campbell and Steenburgh 2017). The main difference between the WRF-LES domain and d03 is in the spatial distribution of the largest precipitation amounts on Tug Hill. There is a banded region where accumulations exceed 60 mm in d03 (Fig. 3b, cyan contour), similar to the simulation by Campbell and Steenburgh (2017). In contrast, this region has a greater north-south extent in the WRF-LES domain, in better agreement with the radar QPE (Fig. 3c).
During the event, the snowband varied in appearance on radar. At 0100 UTC 11 December 2013, snow detected by the KTYX radar began to intensify and a coherent band began to form. By 0220 UTC 11 December 2013, the band appeared rather contorted, and maximum reflectivity values of 30–35 dBZ were observed onshore and on the windward side of Tug Hill (Fig. 4c). Sixteen and a half hours later, at 1850 UTC 11 December 2013, the KTYX radar displayed a more linear primary snowband that extended from Lake Ontario, over the leeward side of Tug Hill, to the Ha-De-Ron-Dah Wilderness (Fig. 4f). Reflectivity values of 30–35 dBZ angled to the northeast across Tug Hill because the shoreline orientation along the southeastern side of the lake (near Oswego, New York) produced land-breeze convergence there in a southwest-to-northeast orientation (see Fig. 7a and Steenburgh and Campbell 2017). A much weaker, second band can be seen to the north of the main band, and this second band formed due to a convergence zone along the northern shoreline of Lake Ontario (see Fig. 7a and Steenburgh and Campbell 2017). Convective cells with reflectivity values of 30 dBZ were embedded in the band at 1850 UTC over the lake and just onshore. A continuous region with reflectivity values of 30 dBZ covered the top of Tug Hill south of the radar location.
The d03 (Figs. 4a,d) and WRF-LES (Figs. 4b,e) reflectivity fields show similar band morphology: by 0220 UTC, the simulated bands are contorted and widen onshore and are slightly north of the observed band, and by 1850 UTC the simulated bands become linear. By 1850 UTC, the band is shifted southward in both domains compared to the observations, which also occurred at this time in the simulation analyzed by Bergmaier et al. (2017). The model-derived reflectivity field is calculated [see Eq. (A1)] to optimally represent the KTYX observations: we assume a 10-cm radar wavelength, and a 0.5° tilt angle from the location and elevation of KTYX (radar range–height circles are shown in Fig. 4e). The same microphysics scheme is used in both d03 and the WRF-LES domain, and reflectivity is calculated using the same equation and the same microphysical inputs. Differences in the reflectivity values when comparing the domains are attributed to differences in the microphysical evolution occurring in each domain.
The main difference between the WRF-LES domain and d03 is that the WRF-LES domain (d04) captures the band’s mesoscale morphology and convective nature better than the 1.33-km simulation (d03). Convective cells with reflectivity values of 30 dBZ are embedded in the band in the WRF-LES domain. In contrast, the d03 band has one continuous region that extends from the lake to the windward side of Tug Hill where reflectivity values are greater than 30 dBZ. The band simulated in d03 is more intense than in the WRF-LES domain with widespread reflectivity values greater than 30 dBZ. The d03 band is similar in appearance to bands simulated by Campbell and Steenburgh (2017, their Fig. 3f) and Bergmaier et al. (2017, their Fig. 2b), which both used the same simulation with an inner domain with 1.33-km horizontal grid spacing.
The coherent, banded region of reflectivity in d03 is a persistent feature. Domain 3 has an extended region from the lake to Tug Hill where reflectivity values are 20 dBZ or greater 70% of the time (Fig. 5a) and a narrow, banded region where reflectivity values are 30 dBZ or greater 30% of the time (Fig. 5d). Frequencies of reflectivity values greater than 20 dBZ from both the KTXY radar (Fig. 5c) and the WRF-LES domain (Fig. 5b) increase going onshore with maximum values corresponding with the highest terrain on Tug Hill. Additionally, there is a broad region in the observations over Tug Hill where 10%–30% of the time reflectivity values are 30 dBZ or greater (Fig. 5f). This region is shifted to the west in the WRF-LES domain (Fig. 5e). The largest reflectivity values more frequently occur over Tug Hill in the observations and the WRF-LES domain, whereas in d03 they more frequently occur along the band and closer to the lake shore.
Probability density plots of reflectivity values at both Sandy Creek (Fig. 6a) and North Redfield (Fig. 6b) confirm that d03 has a high bias in reflectivity values, whereas the probability of reflectivity values from the WRF-LES domain agrees better with the observations, especially at values greater than 25 dBZ. Probability density plots are created using a 13.3-km by 13.3-km square around both locations in both domains, though the results hold using a single grid cell area. An approximate 13.3-km by 13.3-km scan area is used for the radar data. At Sandy Creek, the observed and WRF-LES domain peaks in the probability of reflectivity values are between 20 and 25 dBZ. Additionally, at Sandy Creek, the probability of experiencing a reflectivity value of 32 dBZ is four times higher for d03 compared to observations. This implies that d03 has a mean particle size that tends to be biased high at Sandy Creek. At North Redfield, the probability of experiencing a reflectivity value of 32 dBZ is almost two times higher for d03 compared to observations. Additionally, at North Redfield, the observed and WRF-LES domain peaks in the probability of reflectivity values are between 25 and 30 dBZ, whereas d03 has a peak between 30 and 35 dBZ.
The above analysis shows one of the impacts on the snowband structure that occurs when reducing the model horizontal grid spacing. The horizontal grid spacing used in d03 is too large to resolve the convective cells embedded in the band. Convective cells on radar are approximately 1–2 km in diameter as a best estimate. The convective cells that are resolved in the WRF-LES domain and not in d03 impact the dynamical and ultimately the microphysical evolution of the snowband as seen in the reflectivity field.
The main snowband forms along the southern shore land-breeze front (LBF1; Steenburgh and Campbell 2017) as shown in Fig. 7a. The vertical air motion wair at 1 km AGL is averaged from 1600 to 1800 UTC (during which the band is linear and nearly stationary) for the two domains. The band-averaged updraft is stronger in d03, particularly just onshore (2.75 m s−1, Fig. 7b versus 2.2 m s−1, Fig. 7d), even though updrafts within individual cells are stronger in the WRF-LES domain.
In contrast to what occurs along the main band, the midlake convergence zone (dashed line between LBF1 and the CZ in Fig. 7a) and the northern shore convergence zone (CZ) are weaker in d03 than in the WRF-LES domain. The band of vertical motion over the middle of the lake in the WRF-LES domain is not seen in d03. This midlake band is also seen as a convergence feature at 1800 UTC in the simulations of Steenburgh and Campbell (2017), and it produces an updraft of 0.5 m s−1. This convergence zone is not related to a land-breeze front (Steenburgh and Campbell 2017); it is a center-of-the-lake solenoidal circulation. These convergence zones, the v-wind components, and the averaged vertical motions are stronger in the WRF-LES domain over the lake (Figs. 7a,c). These convergence zones produce narrow updrafts in the WRF-LES domain, and theses updrafts are perhaps too narrow to be resolved on average in d03, especially the center-of-the-lake circulation. The onshore merger of the main band, the northern shore convergence zone (CZ) and the midlake circulation (along with frictional convergence) produces a wide region onshore with updrafts that are stronger in the WRF-LES domain (Figs. 7b,d).
The prominence of a main, dominant band-scale circulation in d03 exists both over the water (leg 3, Fig. 8a) and just onshore (leg 4, Fig. 8c). In fact, the d03 circulation strengthens onshore (Fig. 8c) due to the merger of several convergence zones there (Fig. 7a). On average, the main band is narrower in the WRF-LES domain over the lake compared to in d03 (Figs. 8a,b, compare heating rates and red-shaded regions). In contrast to what happens in d03, the strong, narrow convergence zone over the lake in the WRF-LES domain collapses onshore and merges with the other convergence zones north of it (Fig. 8d). The collapse of convection onshore, consistent with findings from Welsh et al. (2016), is seen in the WRF-LES domain as the largest vertical air motions and heating rates disappear roughly 10 km onshore (Fig. 8d). The collapsed convection from the main band remains upright and deeper than the merged updrafts north of the main band.
The dynamical picture shown above supports the radar analysis: A strong main band dominates in d03, which is conducive for the growth and collection of ice to form large aggregates along the band. In contrast, the stronger, narrower convective elements in the WRF-LES domain break up the band-scale circulation, which appears to hinder the aggregation process as seen by the reduced reflectivity values in that domain (Figs. 4–6). The merger and collapse of convection onshore in the WRF-LES domain helps explain why the concentrated area of high reflectivity values in convective cores expands onshore in areas when the cells collapse (see Fig. 4e).
b. Finescale band evaluation from aircraft observations
The radar evaluation highlights band-scale microphysical differences between d03 and the WRF-LES domain. Aircraft observations and the WRF-LES domain are used to explore the finescale dynamics and how they impact the microphysics.
The strong updrafts in the band produce liquid water which can freeze and subsequently grow into ice particles or remain supercooled and be a source of rime. At 1.7 km MSL over the lake, liquid water contents (LWCs) increase with vertical air motion (wair, Fig. 9a). Both d03 and the WRF-LES domain show this increase, though only the WRF-LES domain captures both the high spread in air motion at LWCs greater than 0.8 g m−3 and the highest observed air motion values (greater than 6 m s−1).
Onshore, both domains have higher LWC values than observed (Fig. 9b). Additionally, the model generally has lower vertical air motion values for the same LWC when LWC values are between 0.2 and 0.4 g m−3. The WRF-LES domain captures the high density and spread of observations with LWCs less than 0.2 g m−3 and vertical air motions less than 2 m s−1.
The hydrometeor vertical velocity field wdBZ (a combination of vertical air motion and reflectivity-weighted particle fall speed) was measured directly by the WCR (W-band). This field is used to analyze the finescale dynamics (and microphysics) of the snowband. To compare the model hydrometeor vertical velocity to the WCR wdBZ, reflectivity-weighted fall speeds (see Molthan et al. 2016) are output for each of the three ice species in ISHMAEL microphysics. A reflectivity-weighted average (using the reflectivity factor, S-band) of these three fall speeds is then computed as a total average value of the reflectivity-weighted fall speed. This value is combined with the vertical air motion wair from WRF to compute wdBZ from the model.
Distributions of wdBZ from the model are sampled along each flight leg (see Fig. 2b) at 1 km AGL from both 1600–1800 UTC (model time), when the simulated band is closer in location to where the observed band was at 1850 UTC, and 1800–2000 UTC (model time), when the simulated band is slightly farther south than the observed band at 1850 UTC (see Fig. 4). During both time periods, the band is linear. Distributions of WCR wdBZ for each flight leg contain transects flown between 1905 and 2029 UTC. Four transects were used for leg 3 (over water), three for leg 4 (over the western foothills of Tug Hill), and one for leg 5 (over Tug Hill).
The WCR wdBZ distribution becomes progressively more narrow going onshore (Figs. 10a–c, black line), and this onshore change in dynamics is corroborated by two ground-based profiling radars located at the eastern shore of the lake and North Redfield (Minder et al. 2015). The peak of this distribution is near −1 m s−1 for each leg which corresponds to ice particles falling at 1 m s−1 (assuming the peak in air motion is at 0 m s−1), the approximate fall speed of aggregates (Locatelli and Hobbs 1974). Over the lake (leg 3) at 1 km AGL, updrafts of over 7.5 m s−1 were measured, and updrafts of 10 m s−1 were measured over the lake up to 3 km MSL (Welsh et al. 2016). The distribution of WCR wdBZ over the lake (Fig. 10a) is positively skewed: downward wdBZ values are not as large in magnitude as upward values. This is characteristic of surface buoyancy-driven boundary layer moist convection (Zhu and Zuidema 2009; Ghate et al. 2010; Lamer and Kollias 2015). Onshore the vertical velocity distribution becomes less skewed and narrower with values ranging from approximately ±3 m s−1 over the windward side of Tug Hill (leg 5, Fig. 10c).
Compared to the WCR wdBZ, the WRF-LES domain captures the distribution of the up and downdrafts from 1600 to 1800 UTC along legs 3 (Fig. 10a, orange line) and 4 (Fig. 10b, orange line), whereas d03 does not capture the larger values of the up and downdrafts over the lake, with peak updrafts of approximately 3 m s−1 (Fig. 10a, blue line), in agreement with Bergmaier et al. (2017). The wdBZ distributions are too narrow for both d03 and the WRF-LES domain compared to observations for leg 5 (Fig. 10c), though the WRF-LES domain has a wider wdBZ distribution in better agreement with observations. The discrepancy between the observations and the WRF-LES domain along leg 5 are likely caused by the complicated nature of the band over land. The combination of the orography (Campbell and Steenburgh 2017), the land-breeze fronts (Steenburgh and Campbell 2017) and the collapse of convection (Welsh et al. 2016) all complicate the onshore dynamics. The WRF-LES domain compares better to the observations at 1600–1800 UTC than at 1800–2000 UTC (Figs. 10c,f) when the simulated band is closer in location (farther north) to the observed band. Thus, the distribution of wdBZ over Tug Hill (leg 5) is likely still influenced by convection (or the collapse thereof) in the band.
The WCR wdBZ spectra are compared to values computed from the WRF-LES domain over both the lake and onshore at 1 km AGL (Fig. 11). The wdBZ spectra are computed from model output using the same wdBZ values computed for Fig. 10. Power spectra are computed along two tracks, corresponding with flight legs 3 (Fig. 11a) and 5 (Fig. 11b), at each model output time from 1600 to 1800 UTC, and then those spectra are averaged and plotted. This is repeated for the time period 1800–2000 UTC. The sampling frequency of the modeled spectra is calculated assuming that an aircraft is flying through the domain at 100 m s−1. Over the lake, the model has similar power spectral densities comparing 1600–1800 and 1800–2000 UTC (Fig. 10c), but over land (leg 5), the power decreases at the later time in agreement with weaker updrafts (Fig. 10f). The inertial subrange is appropriately characterized by the WRF-LES domain down to about 1-km (6–7 Δx, which is the model’s effective resolution; Skamarock 2004) over both the lake and land.
The dynamics-microphysics coupling is further explored by analyzing how reflectivity Z varies with hydrometeor vertical velocity wdBZ. Specifically, this reveals the coupling between updrafts which produce ice and the location of the largest ice particles. WCR Z–wdBZ frequency plots (from 0 to 1 km AGL) along legs 3 (Fig. 12c) and 5 (Fig. 12f) reveal a negative correlation between Z and wdBZ over the lake (leg 3, in the magenta box) and a much weaker but still negative correlation over land. The negative correlation between Z and wdBZ over the lake implies that strong updrafts contain fewer larger ice particles which are being lofted from the tops of these strong updrafts (bounded weak echo regions, BWERs) in a “fountain effect” (see Welsh et al. 2016, their Figs. 7c and 8c) and (see Bergmaier et al. 2017, their Figs. 6a,b). Over land, a much weaker correlation between Z and wdBZ exists. This weak correlation occurs because faster falling particles are generally larger and have higher reflectivity values (not shown, but also seen when air motion is removed and only fall speed is considered in Fig. 12b).
Similar Z–wdBZ frequency plots created from model output over 1600–2000 UTC and 0–1 km AGL to attain a large sample size for both d03 and the WRF-LES domain are shown in Figs. 12a and 12b. The model reflectivity (assumed 10-cm wavelength) cannot be directly compared to the WCR (3-mm) radar, especially since the WCR (a W-band radar) reflectivity starts to plateau around 10–15 dBZ on account of Mie scattering and path-integrated attenuation (Kollias et al. 2007; Matrosov and Battaglia 2009; Welsh et al. 2016). Nevertheless, similar results are expected from the model output: sufficiently strong updrafts should loft ice particles and be associated with lower reflectivity values.
The WRF-LES domain shows a weaker negative Z–wdBZ correlation over the water (Fig. 12b, magenta box) than the observations, but a stronger correlation than d03 (Fig. 12a). Additionally, the WRF-LES domain has a larger spread in wdBZ over both the lake and land compared to d03 (previously shown in Fig. 10). The difference between d03 and the WRF-LES domain is that the average wdBZ decreases in the WRF-LES domain between Z values of 0 to 18 dBZ, in a more similar fashion compared to the observations. This is evidence that the coupling between the microphysics and dynamics in the WRF-LES domain is working in the right direction for Z < 20 dBZ, where the largest of these particles are pushed out of the tops of the strongest updrafts over the lake in a fountain effect. Over land, the updrafts in both d03 and the WRF-LES domain are weaker and there is less of a Z–wdBZ correlation, in agreement with observations. As noted, the wdBZ distribution is too narrow over land for the WRF-LES domain; this is also seen in the Z–wdBZ analysis. The WRF-LES domain and d03 over both the lake and land show a spike in wdBZ corresponding with Z = 25 dBZ. This spike is likely caused by aggregates, which also tend to correspond with positive (upward) wdBZ. Particles still grow and collect in updrafts, but particles are also pushed to the tops of the strongest resolved updrafts in the WRF-LES domain.
In addition to evaluating the simulation based on WCR wdBZ, another way to evaluate the simulated microphysical evolution of the band is through a model-observation comparison of ice particle size distributions. ISHMAEL is a bulk microphysics scheme with three ice species, all of which have fixed gamma distribution shape parameters of ν = 4 (see Jensen et al. 2017). Because the size distribution shape is fixed, we expect differences between the model and the observations. We use this evaluation to confirm that the modeled size distributions evolve appropriately in response to the dynamics and microphysics including sedimentation in the snowband.
The three ice species are combined by binning the size distributions using 200 bins in maximum diameter (D) space. A distribution is created for each model grid cell along a given flight leg from 1600 to 1800 UTC. An ice species must have a mass mixing ratio of greater than 0.001 g kg−1 to be included. For each D bin, the spread between the 25th and 75th percentiles is shaded in Fig. 13.
Over the lake, observations show a general trend of the largest particles (>3 mm) becoming more numerous with decreasing height (Fig. 13a, squares), and both d03 (Fig. 13a, shaded regions) and the WRF-LES (Fig. 13c, shaded regions) domain shows a similar evolution. This suggests that the aggregation process is active over the lake, especially since 1-cm particles are observed there. Both domains have a high bias in number concentrations at D = 1 mm and a low bias in number concentration for particles with D = 0.1 mm compared to observation, in part due to the fixed size distribution shape.
The observations reveal that the number of large particles increases at 1.7 km MSL and just onshore (leg 4) in relatively strong updrafts (wair > 1 m s−1) compared to the downdrafts (wair < 0 m s−1, Fig. 13b). The WRF-LES domain shows a similar result with a distinct separation in ice size distributions between the strong updrafts and all downdrafts (Fig. 13d). There are fewer resolved updrafts in d03 with wair > 1 m s−1 compared to the WRF-LES domain, and the particles in these updrafts in d03 that are larger than D = 5 mm are too numerous compared to observations (Fig. 13b) and could contribute to the high bias in d03 reflectivity in the band. Ultimately, these onshore updrafts contain large particles (aggregates) with relatively low fall speeds which can be deposited downwind and continue to grow by vapor deposition and aggregation over Tug Hill (Campbell and Steenburgh 2017).
The comparison between the model and aircraft observations demonstrates that the WRF-LES domain captures the strongly forced dynamics that occurs in the band better than the 1.33-km domain (d03). Additionally, the WRF-LES domain captures the general evolution of ice particle size distributions and the interaction between ice particles and the dynamics.
4. The impact of better-resolved lake-effect band dynamics on ice particle properties and precipitation
The WRF-LES domain dynamics and microphysics compare reasonably well to radar and aircraft observations. Thus, we explore the ice particle properties including the masses, sizes, shapes, fall speeds, densities, and number concentrations to determine the impact of the in-band resolved convection on the microphysics. It is expected that the lower frequency of high reflectivity values seen in the WRF-LES domain (Figs. 5 and 6) results from the stronger updrafts that produce higher ice number concentrations, more riming and fewer aggregates.
Larger values of ice number concentrations are produced in the WRF-LES domain from the stronger updrafts along both legs 3 and 5 [nucleation, from DeMott et al. (2010), depends on temperature, large aerosol concentration and supersaturation with respect to liquid], particularly for ice-one (nI1), which is planar-nucleated ice, compared to d03 (Table 2). Both domains produce similar ice-one mass mixing ratios (qI1). Ice-two mass and number concentrations are small because ice-two (columnar-nucleated ice) initiation occurs near −7°C, which is near the surface (see Fig. 14) and generally lower in elevation than where most of the ice nucleation occurs for this case. Aggregate (ice-three) mass concentrations are larger in d03 compared to the WRF-LES domain along both legs as expected.
Data from a disdrometer at Sandy Creek were used to derive a median particle diameter throughout the event of 1 mm (Minder et al. 2015). These same data produce a median number concentration of 5.7 L−1. The WRF-LES domain (over 22 h) has a median total ice number concentration of 6.4 L−1, in better agreement with the observations compared to d03 (1.4 L−1). The WRF-LES domain overestimates the median of the total ice number concentration by 12%, but d03 underestimates it by 75%. The resolved dynamics in the WRF-LES domain produce larger ice number concentrations at Sandy Creek in better agreement with observations compared to d03.
A closer examination of the microphysics and dynamics at 1510 along a cross section through Sandy Creek and 1600 UTC along leg 4 (just east of Sandy Creek), reveals that a main cross-band circulation exists in d03 (Figs. 14a,c). Two times and two locations are plotted to show the evolution of the band as it evolves going onshore and in time. The WRF-LES domain has a main circulation in which stronger updrafts are embedded (Figs. 14b,d). The main circulation seen in d03 is conducive to aggregate formation at 1600 UTC (Fig. 14c, high reflectivity values and the magenta line), whereas the main circulation is broken up in the WRF-LES domain and the stronger, more numerous updrafts support larger ice number concentrations through increased nucleation (Figs. 14b,d, white lines), resulting in pockets of aggregates (Fig. 14d, magenta lines) and a more broken reflectivity field. Additionally, the “fountain effect” is evident in the WRF-LES domain at 1600 UTC near y = 21 km and z = 1–1.5 km (Fig. 14d, thick black line). Here the updraft (air motion) is greater than 5 m s−1 and the reflectivity values in the updraft are lower than those above it.
In addition to higher ice number concentrations, larger areas of riming (and more isometric ice particles) occur in the WRF-LES domain (Fig. 15b) compared to d03 (Fig. 15a). The main cross-band circulation that occurs in d03 produces rimed (Fig. 15a, white contour) and comparatively more isometric particles (with ϕI1 = 0.5) near the top of the main updraft and more eccentric (smaller ϕ values in this case) vapor-grown particles elsewhere. These vapor-grown particles in d03 attain maximum diameters DI1 of up to 2 mm near the surface (Fig. 15c), attain fall speeds of less than 1 m s−1 (Fig. 15c, black lines) and have densities of 400–500 kg m−3 (Fig. 15a, black lines) which is typical of vapor-grown, branched, planar particles. Not surprisingly, this is where the aggregate mass mixing ratios and reflectivity values are the largest (Fig. 14a).
In contrast, the narrower, stronger updrafts in the WRF-LES domain cover a larger areal extent across the band than in d03 at 1510 UTC. This supports higher ice number concentrations aloft (Fig. 14b) and more regions where riming can occur (Fig. 15b, white contours). These rimed particles have ϕI1 = 0.5–1.0 (Fig. 15b), densities of 300–400 kg m−3 (Fig. 15b, black contours), DI1 = 0.5–1 mm (Fig. 15d) and ice particle fall speeds greater than 1 m s−1 (Fig. 15d, black contours) at some locations. Rimed particles were observed at the tops of the strongest updrafts at 3 km MSL (Welsh et al. 2016); riming occurs in the WRF-LES domain at the tops of the strongest updrafts and to a much lesser extent in d03.
The impact of the resolved convective cells on the microphysics are twofold. First, convective cells hinder the aggregation process by increasing both ice number concentrations through increased nucleation rates and riming rates, thus limiting increases in mean particle size and reducing reflectivity values along the band. Second, convective cells include narrow, strong up and downdrafts. These convective cells likely disrupt the band-scale circulation (see Bergmaier et al. 2017); the WRF-LES domain band is thus composed of convective cells in a banded orientation.
Nucleation generally occurs near the top of the band (see Fig. 14), and produces, small, spherical ice in the model which reduces the bulk size of an ice distribution and makes the ice particles more spherical on average. In ISHMAEL, a two-moment bulk microphysics scheme, nucleation (adding small particles) can reduce the mean size of a distribution. Physically, adding particles to a distribution would not impact in situ aggregation rates. Regardless, these newly nucleated particles will need time to grow by vapor deposition to sizes that can collect.
Riming causes ice particles to become more spherical. Smaller, more spherical particles collect with a smaller efficiency than larger, more eccentric particles (Connolly et al. 2012). Nucleation and riming shift the distribution of reflectivity values toward lower ones in better agreement with observations (Fig. 6). The orographic and stratiform lift over Tug Hill is a weaker forcing than convective cells, which provides an environment with significant vapor growth, less riming and more aggregation (Campbell and Steenburgh 2017). Thus we expect and do see the highest reflectivity values over the terrain in the WRF-LES domain.
At Sandy Creek for the duration of the event, a smaller percentage of aggregates existed at the lowest model level in the WRF-LES domain along with a higher percentage of stronger updrafts than in d03 (Table 3). In the WRF-LES domain, 2.8% of the lowest model level ice at Sandy Creek is ice-two and this ice is all quasi spherical (0.8 < ϕI2 < 1.2). Seven percent of the lowest model level ice (considering total ice) at Sandy Creek in the WRF-LES domain is ice-one with ϕI1 > 0.25 and 1% of all lowest model level is ice-one with ϕI1 > 0.5. The averaged mass-weighted diameters and mass-weighted fall speeds of this ice, partitioned by ϕ, are given in Table 3. Using data from Locatelli and Hobbs (1974), densely rimed dendrites with maximum diameters of 1.7 mm fall at 0.7 m s−1, though note that the range of their data starts at a diameter of 1.8 mm. Graupel particles with densities of 200–450 kg m−3 and maximum diameters of 0.7 mm fall at 1.3 m s−1. Thus modeled ice particles with ϕI1 > 0.1 have average maximum diameter and fall speed values that match those of observed densely rimed particles. The average maximum diameters and fall speeds of particles with ϕI1 > 0.5 match those of observed graupel. Modeled particles with ϕI1 > 0.25 have average fall speeds that are between the observed fall speeds of densely rimed particles and graupel. Based on this analysis 3.8% of ice at the lowest model level at Sandy Creek in the WRF-LES domain is graupel and 9.8% is densely rimed (including graupel) based on particle properties. This agrees well with results from Welsh et al. (2016), who used disdrometer data to estimate that 10% of the precipitation at Sandy Creek was graupel-like based on fall speeds, and the most probable particle size was a diameter of 0.5 mm.
Ultimately, the pockets of rimed particles can produce narrow but significant increases in the vertical flux of ice near the surface at Sandy Creek (Fig. 16). While rimed particles may only account for 10% of ice at the surface (for this case), those particles may help significantly increase precipitation rates locally for brief periods of time.
The majority of riming over 22 h in d03 (Fig. 17a) is associated with the main circulation (LBF2, see Fig. 7). Riming is associated with LBF2 for the WRF-LES domain (Fig. 17b), but also with the resolved convective cells over the lake. There is a larger area of riming over the lake in the WRF-LES domain compared to d03. The additional offshore riming in the WRF-LES domain results in a lower flux of aggregates just onshore on average in the WRF-LES domain compared to d03 (Figs. 17a,b, blue contours).
The 5-min gauge data (Steenburgh et al. 2014b) shows that precipitation rates exceeded 4.23 mm h−1 (2 in. h−1 of snowfall assuming a 12–1 snow-to-liquid ratio, a good estimate at Sandy Creek using data from 0600 to 1200 UTC) at Sandy Creek for 145 min. At Sandy Creek and over the lake (west of Sandy Creek), liquid equivalent precipitation rates estimated from Level II radar QPE exceed 4.23 mm h−1 for 30–60 min (Fig. 18c). The frequency of these heavy precipitation rates, which represent the most hazardous part of lake-effect snow storms, are matched better in the WRF-LES (60–90 min, Fig. 18b) compared to the gauge data than in the d03 simulation (30–60 min). Also, they are seen in the same location in the WRF-LES domain (Fig. 18b) compared to the radar QPE but are not seen in d03 (Fig. 18a), where precipitation rates this large and for at least 30 min do occur just onshore. At North Redfield, precipitation rates exceeded 4.23 mm h−1 for 220 min according to the gauge data (Steenburgh et al. 2014a), and this amount is in slightly better agreement with the WRF-LES domain (210–240 min) than d03 (240–270 min). There is a large, banded region in d03 where precipitation rates exceed 4.23 mm h−1 for 6.5–7 h, which is not seen in the radar observations. In d03, the duration of relatively heavy snowfall is missed just onshore and is over predicted for a banded region east of North Redfield, when verified against the radar QPE.
The OWLeS IOP2b case is simulated using a nested WRF configuration with the innermost domain utilizing 148 m horizontal grid spacing, which is ninefold smaller than used in previous simulations of the case (Bergmaier et al. 2017; Campbell and Steenburgh 2017; Steenburgh and Campbell 2017). Results using this high-resolution WRF-LES domain are compared with those using a coarser 1.33-km horizontal grid spacing domain (d03). A direct result of the increased resolution is that the model is able to capture the strongest updrafts in the band, up to 7.5 m s−1, in better agreement with aircraft observations of Doppler hydrometeor vertical velocity.
There are several changes that occur in the microphysical evolution of the band when stronger updrafts are resolved in the WRF-LES domain. These changes are apparent when comparing the simulations to radar, surface and aircraft observations. Stronger updrafts lead to increased ice nucleation rates and riming rates. These increases in ice number concentrations and rimed particles produce a higher fraction of ice that is not purely vapor grown, and because these more spherical particles collect with a lower efficiency than vapor grown ones, aggregation rates are reduced. Radar observations generally support this microphysical picture of the band at two locations onshore. The most probable reflectivity value onshore is less in the WRF-LES domain than the 1.33-km domain because aggregation rates are reduced, and because smaller, more numerous particles populate the band. Aircraft and surface observations corroborate this idea. Additionally, the stronger updrafts disrupt the cross band circulation, which is on average weaker onshore in the WRF-LES domain.
The overall realism of the d03 simulation in terms of its ability to capture the precipitation as well as the general reflectivity field implies that the band dynamics and the evolution of the microphysics including the impact on the reflectivity field are dominated by the band-scale (mesoscale) convergence (Bryan et al. 2003). Nonetheless, differences between d03 and the WRF-LES domain in the reflectivity field highlight how better resolving lake-effect dynamics impacts the microphysics, including how stronger updrafts change the microphysics. Understanding the impact of better resolved dynamics on microphysics is important because both operational and research models are being run at increasingly higher resolution. The WRF-LES domain can be used also as an assessment of radar-based precipitation estimation in lake-effect snow storms. Radar-derived quantitative precipitation estimates in general tend to underestimate heavy precipitation in wintertime over land because of radar beam overshoot (Cocks et al. 2016). Additionally, misclassification of precipitation type, faulty Z–S relationships, and sublimation/growth of ice below the radar beam can all lead to errors in radar-derived precipitation estimates. Continued high-resolution modeling studies and targeted field campaigns can better constrain these precipitation estimates, which are critical for real-time warnings.
As models are run at increasingly higher resolutions, a need to change or improve model parameterizations, such as microphysics, may be necessary. As found in this study, a shift in microphysical process rates occurs when more intense updrafts are resolved. It is beyond the scope of this work to determine how traditional microphysics schemes, which use predefined, discrete categories such as cloud ice, snow, and graupel, will handle higher resolutions for modeling lake-effect snowbands. Traditional schemes must formulate conversion processes between categories, and these processes are often ad hoc and based on thresholds that can lead to large, discrete changes in simulated clouds and precipitation (Morrison and Grabowski 2008). Thus, we hypothesize that traditional microphysics schemes may be more sensitive to changes in resolution; better resolving updrafts will increase riming rates which could tip snowband simulations from being snow-dominated to graupel-dominated. Schemes like ISHMAEL microphysics and P3 (Morrison and Milbrandt 2015) that eschew traditional ice categories and smoothly evolve ice particle properties like density and fall speed may be more adept to handle increased model resolution. Testing this hypothesis is left to future work.
The authors thank the National Science Foundation’s NSF-AGS postdoctoral research fellowship program (Grant AGS-1524267), which was instrumental for this research. B. Geerts and P. Bergmaier acknowledge NSF Grant AGS-1258856 and P. Bergmaier acknowledges NSF Grant DUE-1821566. Implementation of the high-resolution terrain was made possible by the one and only Michael G. Duda. Discussions with Peter Sullivan about turbulence were enlightening. We are thankful for useful comments by James Pinto. We thank the UWKA flight crew for their hard work. We thank everyone collecting data during the OWLeS campaign. NCEP Reanalysis data were provided by NOAA/OAR/ESRL PSD, Boulder, Colorado, from their web site at https://www.esrl.noaa.gov/psd/. The National Center for Atmospheric Research is sponsored by the National Science Foundation. The authors are thankful for comments from three anonymous reviewers.
Jensen ISHMAEL Microphysics Scheme Details
a. Scheme overview
All three ice species in ISHMAEL are modeled using spheroids, defined by an aspect ratio ϕi = ci/ai, where ci and ai are the spheroidal semiaxis lengths. A spheroid is oblate (prolate) if ai > ci (ci > ai). For planar- and columnar-nucleated ice, four prognostic variables are predicted: mass, number, volume, and volume times aspect ratio mixing ratios (two volume-based mixing ratios are needed to predict two axis lengths). Aggregate shape and density are currently assumed (see Jensen et al. 2017), and therefore, only mass and number mixing ratios are predicted for aggregates. Microphysical processes such as vapor growth and riming cause the bulk shape (habit) and density of planar- and columnar-nucleated ice to evolve. Ice particles thicken (become more isometric) when they collect rime, and this thickening is parameterized in the scheme through a change in the bulk shape of ice which additionally causes the fall speed to increase (Jensen and Harrington 2015). Each ice species uses a fixed gamma distribution shape parameter of ν = 4.
The spheroidal axes are related by a power law (Harrington et al. 2013). The parameter is a time-integrated quantity and tracks the shape history, while is chosen such that ice is spherical at nucleation (Harrington et al. 2013). The mass of a spheroidal ice particle is , where ρi is the effective density and accounts for air in the ice (e.g., space between dendritic branches, hollow areas in columns and air pockets between rime). Using the power-law relationship between ci and ai, mass can be written as , where and . Thus, mass can be written using the traditional mass-dimensional formulation, but βi changes as ice particle shape changes in ISHMAEL. In traditional bulk microphysics schemes βi is fixed based on ice category (e.g., βi = 3 for graupel).
b. Reflectivity equation
The equivalent reflectivity factor Ze is calculated assuming a 10-cm radar wavelength and Rayleigh scattering only (the radar cross section becomes a function of mass squared) and accounts for the reduced reflective capacity of ice compared to liquid (Smith 1984; Ferrier 1994; Koch et al. 2005; Hogan et al. 2006; Westbrook et al. 2006; Bryan and Morrison 2012). For a single ice species,
where Ki = 0.176 is the dielectric constant for ice, ρs is the density of solid ice, and the mass mi depends on the shape and density evolution as shown above.
c. Number-weighted aspect ratio equation
This article is included in the Ontario Winter Lake-effect Systems (OWLeS) Special Collection.
The Jensen ISHMAEL microphysics is publicly available as of WRF, version 4.1, mp_physics = 55.