*2.1. Composite Analysis*

Composite surface and upper-air patterns were generated using the NARR [24]. Though on a native 32 km horizontal grid, this dataset was averaged to a lower resolution, 16 × 16, and 1.25◦ (longitude) by 0.94◦ (latitude) grid centered on the Grand Forks NWSFO CWA. This was done to reduce the computational cost of the SOM and to facilitate future comparisons to other datasets (e.g., weather or climate model output). While a number of re-analyses are now available, NARR was chosen due to the authors' familiarity with this dataset, along with prior studies that demonstrated favorable performance over the region [25–28]. Given the variables and resolution used, it is anticipated that similar results would be found if other current generation re-analyses were used (e.g., ERA-Interim [29]). This assertion may not be valid in regions with more limited surface and upper-air observations where re-analyses are more poorly constrained.

Using storm data, available surface observations, and the NARR, midpoint times were estimated for each blizzard event. Patterns were composited for the four primary patterns using blizzard midpoint times and +/−12-hr before and after these points. For patterns that contained a mid-latitude cyclone, the minimum mean sea level pressure (MSLP) within the domain was identified and tracked over this time.

#### *2.2. Objective Classification Using a SOM*

To objectively classify atmospheric patterns, the Self-Organizing Map (SOM) [30] technique was used. A competitive neural network, SOMs are most similar to a K-means clustering algorithm [31]. Unlike K-means clustering, SOMs include a neighborhood function during the training process. The result is a topological (feature) map that allows clusters (nodes) to a) span the data space and b) relate to each other in a two-dimensional matrix. This latter property allows users to be less concerned with the exact number of clusters to choose and instead to focus on clusters that are relevant for their analysis purposes. While this alone makes it a useful algorithm for pattern recognition, SOMs hold other advantages (handling of noise, no a priori assumptions of data, better identification of pattern mixing) over common techniques such as Principal Component Analysis (PCA)/Empirical Orthogonal Functions (EOFs) [32–34]. As a result, SOMs are now commonly used in the fields of meteorology and oceanography. For additional information, the reader is referred to earlier surveys of SOM studies [35,36].

The process of creating a SOM follows the strategy employed in earlier work by the author [37], and the reader is referred to this study for more details on the nuances of SOM creation. To summarize the process, a user must first select data for input then reduce the multi-dimensional meteorological data into input vectors that the SOM performs the clustering on. SOMs are trained in a two-step process that first determines the orientation of the feature map then iterates to a final solution that seeks to minimize the error between the training dataset and the final classification of nodes [31]. These stages require the selection of user parameters such as the map size, training length, learning rate, and the neighborhood radius. After the SOM is created, training samples are compared to each node within the feature map and classified to the node with the minimum Euclidean distance.

Consistent with the generation of composite patterns in the previous section, the spatially averaged 16 × 16, 1.25◦ (longitude) by 0.94◦ (latitude) NARR was used to train the SOM. Based on the results of the compositing process, variables that showed significant variability across patterns were used, and these included 500 hPa geopotential heights, MSLP, and surface temperatures. Other combinations of variables were tried, but the inclusion of 500 hPa geopotential heights made the largest difference in the ability of the SOM to segregate patterns. Identical to [37], variables were computed as anomalies from the field mean at each given time step. This allowed the SOM to focus on the gradients in variables, minimizing the issues of biases or variability that vary by season or exist when patterns are compared across multiple datasets (e.g., NARR vs. climate model data) which is useful for future studies. To capture the progression of systems across the domain, each training sample included time steps at the midpoint and +/−12 hrs. With three total variables, a 16 × 16 region, and three times for each case, input vectors used to train the SOM had a length of 2304 elements. All variables were normalized to a common scale to contribute equally to the SOMs. In total, 93 blizzard cases (a subset of 37 winters from 1979–2018 and 2015–2016) were used as input vectors due to the availability of NARR data at the time the SOM was created and as a goal to classify future patterns. Errors (Euclidean distances) for classified patterns in the latter two seasons were similar to the trained data. This suggested that (1) training the SOM with all 100 patterns would not significantly alter the results, and (2) ample variability was captured in the SOM and the methodology is useful for pattern recognition purposes.

A key parameter of the SOM (and other objective classification techniques) is the number of classes chosen. For classifications of atmospheric states, this decision will be dependent on the purpose of the study as well as the number of samples being used to train the SOM. Too few classes can smooth out the details of patterns, while too many will lead to situations where some SOM nodes have no observed patterns classified to them [37]. With a relatively small number of cases (~90–100) and the purpose of comparing the SOM to subjectively classified classes, a rectangular 8-class (4 × 2) is presented. Larger maps were also created, but they did not provide further insight to the results shown herein.

The SOM was generated using SOM\_PAK software, which is freely available [30]. Within this package, the 'vfind' program was used; this program randomly initializes a SOM feature map a specified number of times and selects the map that minimizes the lowest quantization error. Following the guidelines of SOM\_PAK [30], settings for 'vfind' included a training length that increased and learning rates and neighborhood radii that decreased between the two steps in the training process (Table 1).


**Table 1.** Self-Organizing Map (SOM) settings used with the SOM\_PAK command 'vfind'.

#### **3. Results and Discussion**

#### *3.1. General Characteristics*

During the 39-year period, 100 total blizzards were reported, averaging 2.6 events per year. An annual and seasonal breakdown of these events is provided in Figure 3. RRV Blizzards are highly variable, with seasons varying from 0–10 events (Figure 3a). Record years (10 events) included the infamous 1996–1997 winter that concluded with the catastrophic RRV flood [38,39] and the 2013–2014 winter that did not have significant flooding. On the other end of the spectrum, three seasons (1986–1987, 1990–1991, and 2011–2012) did not have any recorded blizzards. While the source of this

variability is beyond the scope of this study, snowfall and cyclone variability in this region have been tied in part to phases of the El Niño Southern Oscillation (ENSO) and the North American Oscillation (NAO) [40,41].

**Figure 3.** (**a**) Annual and (**b**) bimonthly number of *Storm Data* blizzards for the winter seasons of 1979–1980 and 2017–2018. Named blizzards by the Grand Forks Herald are provided by the red dots in panel (**a**).

Out of curiosity, named blizzards from the Grand Forks Herald newspaper were also compared for annual totals. Over the shorter period (1989–1990 and 2017–2018), the newspaper named 63 (vs. 76) blizzards, and the datasets had a correlation of 0.77. Provided that the distribution area for the paper is smaller than the CWA, these results are expected. While some specific years had more events recorded by the paper vs. *Storm Data*, this is attributed to events that were stronger winter storms but did not meet official blizzard criteria.

Blizzards have been reported from October–April, with the bulk of the events occurring from the 2nd half of December to the 1st half of March (Figure 3b). The most frequent period of occurrence was the 2nd half of December, with a total of 17 blizzards over the 39-year period. A unique aspect of the seasonal cycle is the bimodal distribution with a well-defined lull in late February. This is in agreemen<sup>t</sup> with North American cyclone climatologies that indicate a relative minima of cyclones during February [42].

## *3.2. Composite Analysis*

Classifications of the 100 classified blizzards were used to generate composite patterns from the NARR. Of the 100 patterns, two patterns were sufficiently different that they did not fit any of the four categories, and these were omitted from the composite analysis. These included two events driven by southerly winds well ahead of weaker mid-latitude cyclones on 6 March 2014 and 31 December 1996. The remaining composite patterns are now described.

Of the four patterns, Colorado Low blizzards feature the strongest mid-latitude cyclone (Figures 4a–7a) and resemble prior composites of this type [43]. Tracking from Northeast Colorado to North Wisconsin, the composite minimum MSLP decreases from 1002–1000 hPa from 12 h prior to the midpoint of the event. With a storm track south of the RRV, the region is predominately under northerly surface winds that strengthen and shift direction from the East-Northeast to the Northwest as the cyclone progresses eastward. While not shown (and noted earlier), these systems are responsible for the highest snowfall totals, as the region falls within the precipitation shield north of the cyclone track [44]. Aloft, these events are associated with the progression of a well-defined trough that deepens over the region (Figures 5 and 7a). In response to the passage of this trough, strong 500 hPa height falls are found leading up to the event, with the maximum decrease located just northeast of the surface low. In some cases, these troughs are associated with an upper-level closed low, although this definition has been lost to some extent during the compositing process.

**Figure 4.** North American Regional Re-analysis (NARR) composite plots of mean sea level pressure (MSLP) (hPa), surface wind barbs (kts), and surface temperatures (◦C) 12 hr prior to the midpoint of (**a**) Colorado Low, (**b**) Alberta Clipper, (**c**) Hybrid, and (**d**) Arctic Front blizzards. 12-h MSLP change (midpoint—12 hr prior) is provided by shaded contours while composite mean cyclone tracks are denoted by the thick black lines for select classes.

Blizzards associated with Alberta Clippers also feature a well-defined, albeit weaker (1008–1006 hPa) mid-latitude cyclone (Figures 4b–7b). Consistent with the name and prior composites [19], these systems track east-southeast from southern Canada across the RRV with the cyclone center eventually reaching Northeast Minnesota and North Wisconsin. One should note that the composite cyclone tracks appear to be shifted east compared to the Colorado Low and barely encompass the passage of the low out of Canada. Because these tracks were identified symmetrically around the midpoint of the blizzard conditions, this implies that poor visibility primarily occurs after the passage and development of the surface cyclone (Figure 4b). Aloft, conditions leading up to the event feature stronger west-northwest 500 hPa flow with maximum height falls located over Minnesota, just ahead of a developing short-wave trough (Figure 5b). By the midpoint of the blizzard, this trough has amplified and progressed eastward across the domain with 500 hPa winds shifting to the Northwest.

As noted earlier, the Grand Forks NWS defines Hybrid events as those with characteristics of multiple patterns, and this is also true of the composite patterns (Figures 4c–7c). At the surface, this class manifests itself as a mid-latitude cyclone track that begins farther north (south) of a Colorado Low (Alberta Clipper). The minimum pressure and intensity of the wind field are similar to that of the Alberta Clipper, but it has weaker 12-hr pressure rises/falls (Figure 6c). This latter property can be attributed to the slower progression of Hybrids vs. Alberta Clippers. At 500 hPa, Hybrids are associated with weaker, near-zonal flow 12 h prior to the midpoint of blizzard conditions (Figure 4c). Compared to the Alberta Clippers, the short-wave trough is in a similar position, but the orientation of the flow leads to a more neutral tilt. Unlike the aforementioned pattern, Hybrids feature more deepening of the upper-level low/trough by the midpoint of the blizzard, similar to what is seen for Colorado Lows.

**Figure 5.** NARR composite plots of 500 hPa Geopotential heights (m) and wind barbs (kts) 12 hr prior to the midpoint of (**a**) Colorado Low, (**b**) Alberta Clipper, (**c**) Hybrid, and (**d**) Arctic Front blizzards. 12-h 500 hPa height change (midpoint—12 hr prior) is provided by shaded contours.

Arctic Fronts are the final, and arguably most unique, composite pattern identified with RRV blizzards (Figures 4d–7d). Unlike the other patterns, no centralized region of low-pressure is seen at the surface. Instead, this pattern features an elongated southwest to northeast oriented surface trough associated with the developing Arctic Front (Figure 4d). As the event progresses, surface pressures rapidly rise and northerly winds strengthen behind the front, as the Arctic High develops to the Northwest, increasing the gradient in MSLP and Cold Air Advection (CAA). Because the surface trough/Arctic Front is often associated with a more distant cyclone, 500 hPa patterns are more dissimilar from the other patterns (Figures 5 and 7d). These events are characterized by strong northwest flow that eventually develops a southwest trough (passing vorticity maxima) in the eastern half of the domain by the mid-point of the event. As a result, the region ends up residing under large (60 m) height rises associated with a strengthening jet stream and implied Negative Vorticity Advection (NVA). Compared to the Alberta Clippers and Hybrid events, 500 hPa winds associated with Arctic Fronts are approximately double in magnitude (80 vs. 40 kts). This leads to a vertical wind profile (not shown) that implies downward transfer of momentum in a regime of subsidence is a key mechanism for reaching blizzard criteria for winds. The presence of CAA, NVA, and subsidence matches many of the checklist items for impactful post-cold frontal winds [21].

Meteorological patterns responsible for blizzard events have preferred periods of occurrence (Figure 8). Early and late season events (October, November, and April) are primarily due to Hybrid and Colorado Lows, with only one (Alberta Clipper) event not fitting these categories. These classes have bimodal distributions with Colorado Lows (Hybrids) peaking in December and March (January and March), respectively. Alberta Clippers occur from December to March, with the majority of the events occurring during January and February, consistent with [19]. Arctic Fronts, commonly responsible for ground blizzards, are more common during the late winter with events between December to March and a maximum in January. As a result of these distributions, January ends up being the most diverse month with relatively constant fractions (0.2–0.3) across the categories (Figure 8b). As noted earlier, the lull in February is consistent with extra-tropical cyclone climatologies, and this is seen in Figure 8 as a reduction of Colorado and Hybrid lows in this month.

**Figure 6.** As in Figure 4, except for the midpoint of the blizzard. 12 hr MSLP change (12 hr post—midpoint) is provided by shaded contours.

**Figure 7.** As in Figure 5, except for the midpoint of the blizzard. 12-h 500 hPa height change (12 hr post—midpoint) is provided by shaded contours.

**Figure 8.** (**a**) Number and (**b**) fraction of monthly blizzards for the winter seasons of 1979–1980 and 2017–2018, separated by type.

#### *3.3. Objective Classification of Patterns Using a SOM*

Surface and 500 hPa analyses for the midpoint of blizzard events in the eight-class SOM are shown in Figures 9 and 10. The SOM shows a progression of patterns that shift from cold fronts associated with CAA (nodes 1/5) to deeper surface lows (nodes 3/4/7/8). These patterns have 500 hPa analyses similar to those seen in earlier composites. For example, nodes 1/5 resemble Arctic Front patterns with a strong northwesterly flow aloft, while the rightmost nodes appear as Colorado Lows with either upper-level toughing (nodes 3/4) or a closed low (node 7/8). The progression of systems is also similar to the composites shown earlier (not shown). For example, the rightmost nodes (4/8) progress northeastward like a Colorado Low. Shifting from right-to-left, the mid-latitude cyclones become weaker and have tracks that are displaced northerly, consistent with Hybrid/Alberta Clipper type systems. While the SOM has many positive traits when compared to the composites, it by no means is a perfect reproduction of the subjectively classified classes. For example, Arctic Fronts have pressures that are too low 12-hrs prior to the midpoint of events resulting in weak cyclones vs. the open trough seen in Figure 4. This is undoubtedly a result of the neighborhood function within the SOM smoothing this category with other mid-latitude cyclone nodes. Despite this issue, quantitative comparisons of blizzards events to both the composite patterns, and the eight-class SOM yielded better agreemen<sup>t</sup> for the SOM with mean Euclidean distance ~30% lower (315.6 vs. 444.1). Increasing the SOM to a 5 × 3 (or larger) map mitigates this issue and further lowers the mean Euclidean distance; however, this happens at the expense of decreasing the number of blizzards that occur per class (not shown).

**Figure 9.** MSLP (hPa, dashed lines) and surface temperature (◦C, filled contours) anomalies during the midpoint of blizzards for the eight-class (2 × 4) SOM. Nodes are identified by the external numbers ranging from 1–4 (5–8) for the top (bottom) rows.

**Figure 10.** As in Figure 9, except for 500 hPa height anomalies (shaded and dashed contours).

As a final test of the SOM's ability to segregate patterns, the subjectively classified events were categorized into the eight-class SOM (Figure 11). As expected by the meteorological interpretation of the nodes, patterns have distinct areas of occurrence. Colorado Lows (Figure 11a) only occur on the right-hand-side of the SOM, with the majority of the cases occurring within Node 4. Alberta Clippers occur within the left-most six nodes, with most occurring within Nodes 1, 2, 6, and 7. Hybrids, which are subjectively defined as patterns with features of multiple patterns are the only category to occur within every node. That said, the majority of these cases occur within Nodes 3 and 7, in-between Colorado Lows and Alberta Clippers. Finally, Arctic Fronts are primarily on the left-hand-side of the SOM with the majority of the cases occurring in Node 5. From a probability stand-point, categories within the SOM can be arranged in a column fashion, with probability of occurrence shifting from Arctic Fronts (Nodes 1/5), to Alberta Clippers (Nodes 2/6), to Hybrids (Nodes 3/7), to Colorado Lows (Nodes 4/8). By doing this, the time period of occurrence for these categories gives results similar to the results shown in Figure 8 with seasonal occurrence of SOM nodes varying by column (Table 2).

**Figure 11.** Percent of (**a**) Colorado Low, (**b**) Alberta Clipper, (**c**) Hybrid, and (**d**) Arctic Front blizzards identified within each of the eight SOM nodes.



#### *3.4. Discussion and Future Work*

The good agreemen<sup>t</sup> between subjectively and objectively identified blizzard patterns provides evidence that the characteristics of these events, including types of patterns and time periods of occurrence, are well understood. The use of a relatively small SOM and inclusion of only several

variables to obtain this finding is a positive result that suggests SOMs can be used to investigate a number of outstanding questions regarding blizzards. Some of these activities are now discussed.

Within the realm of weather prediction, a constant struggle is determining whether visibility criteria will be met to justify and verify products such as NWS blizzard warnings. Though blowing snow parametrizations exist [45,46], they are not currently included in operational Numerical NWP models within the US. Instead, local forecasters must use empirical models that determine a probability of blowing snow given conditions such as wind speed, air temperature, and snowpack conditions [47]. The context of how these events fit within the scope of forcing mechanism (type of event) is currently only considered subjectively (e.g., Arctic Fronts are harder to forecast than Colorado Lows). A possible solution is for SOMs to provide real-time identification and classification of forecast atmospheric patterns from deterministic or ensemble NWP systems. While this could be done subjectively, burden is placed on the forecaster to identify patterns within the large range of modeling systems now available. Further, there is always the issue of human bias. A hypothetical system can identify the fractional number of ensemble members with forecasted blizzard patterns. Pattern typing can then inform forecasters on the scope of impacts. For a system such as the Global Ensemble Forecast System (GEFS), doing this subjectively would require the forecaster to individually inspect 21 members, a process that is too arduous for an operational setting.

Prior to the implementation of such a system, the SOM must consider null cases to understand the nuances between patterns that do and do not produce blizzard conditions. In practice, this can be done in two di fferent ways. As presented here, a small SOM developed solely from blizzard events can be used by defining a threshold Euclidean distance that counts a pattern as a hit. A caveat with the presented SOM is distances that are still quite large; this would lead to a higher risk of false alarms. Instead, a larger SOM could be used to reduce the threshold Euclidean distance needed. Alternatively, SOMs can be produced using all available patterns during the winter. To encompass the increased variability in patterns (blizzard and null cases alike), these SOMs must be larger. Instead of using a threshold to define events, observed blizzard events can be mapped to the larger SOM to understand which nodes are associated with these events.

Regardless of the exact methods, an interesting avenue of future work is a retrospective analysis on all historical patterns. This can provide insight into events that may have been missed by the observation system or were too limited in scope to fit within the traditional zone/county verification process at the NWS. Pattern recognition can also be extended farther back in time using datasets such as the 20th Century Re-analysis [48] to yield a long-term climatology of blizzards.

How the frequency and intensity of RRV blizzards may change in a warming climate is also unknown. Previous studies have focused on how precipitation or cyclone frequency may change independently. From the Clausius–Clapeyron relationship, a warmer climate will dictate higher amounts of column water vapor and, thus, precipitation [49]. During the winter, however, there will be a balance between warmer temperatures, column water vapor, and precipitation phase. Overall, a general decline in snow cover has been found for the northern hemisphere, and much of this is due to a significant shortening of the snowy season [50,51]. Despite this trend, the RRV region has seen an increase in snowfall, especially for higher end events with 2+ inches [52,53].

Regarding forcing mechanisms for RRV blizzards, mixed results have been found for extratropical cyclones. While some studies sugges<sup>t</sup> a decrease in Northern Hemisphere wintertime cyclone frequency [54,55], other work suggests the strongest cyclones have intensified or could further intensify in future climate projections [56–59]. Concerning specific patterns identified within the present study, there is a projected decrease (increase) in Alberta Clippers (Colorado Lows) over North America [59]. It is unknown how Hybrid lows or Arctic Fronts may change, and this is an avenue of work into which SOMs can provide insight, as they can provide information on type and frequency of occurrence of patterns.
