# Coherent Bimodal Events in Ensemble Forecasts of 2-m Temperature

CAMERON BERTOSSA,ª,ᵇ PETER HITCHCOCK,ᵇ ARTHUR DEGAETANO,ᵇ AND RIWAL PLOUGONVENᶜ

(Manuscript received 20 March 2023, in final form 24 July 2023, accepted 26 July 2023)

**ABSTRACT:** A previous study has shown that a large portion of subseasonal-to-seasonal European Centre for Medium-Range Weather Forecasts (ECMWF) ensemble forecasts for 2-m temperature exhibit properties of univariate bimodality, in some locations occurring in over 30% of forecasts. This study introduces a novel methodology to identify “bimodal events,” meteorological events that trigger the development of spatially and temporally correlated bimodality in forecasts. Understanding such events not only provides insight into the dynamics of the meteorological phenomena causing bimodal events, but also indicates when Gaussian interpretations of forecasts are detrimental. The methodology that is developed allows one to systematically characterize the spatial and temporal scales of the derived bimodal events, and thus uncover the flow states that lead to them. Three distinct regions that exhibit high occurrence rates of bimodality are studied: one in South America, one in the Southern Ocean, and one in the North Atlantic. It is found that bimodal events in each region appear to be triggered by synoptic processes interacting with geographically specific processes: in South America, bimodality is often related to Andes blocking events; in the Southern Ocean, bimodality is often related to an atmospheric Rossby wave interacting with sea ice; and in the North Atlantic, bimodality is often connected to the displacement of a persistent subtropical high. This common pattern of large-scale circulation anomalies interacting with local boundary conditions suggests that any deeper dynamical understanding of these events should incorporate such interactions.

**SIGNIFICANCE STATEMENT:** Repeatedly running weather forecasts with slightly different initial conditions provides some information on the confidence of a forecast. Occasionally, these sets of forecasts spread into two distinct groups or modes, making the “typical” interpretation of confidence inappropriate. What leads to such a behavior has yet to be fully understood. This study contributes to our understanding of this process by presenting a methodology that identifies coherent bimodal events in forecasts of near-surface air temperature. Applying this methodology to a database of such forecasts reveals several key dynamical features that can lead to bimodal events. Exploring and understanding these features is crucial for saving forecasters’ resources, creating more skillful forecasts for the public, and improving our understanding of the weather.

**KEYWORDS:** Nonlinear dynamics; Atmosphere-ocean interaction; Statistical techniques; Ensembles; Probability forecasts/models/distribution; Atmosphere-land interaction

## 1. Introduction

The use of ensemble forecasts in order to create a probabilistic viewpoint of the future state of the atmosphere has been an important development in numerical weather prediction (Toth et al. 2003; Leutbecher and Palmer 2008; Vannitsem et al. 2018). Often, these forecast distributions are treated inherently as Gaussian, where the mean of the ensemble is assumed to be the most likely future state and the variance is the uncertainty associated with this prediction. This assumption, however, is not always appropriate. The present work extends a previous study that demonstrated that bimodality is in fact quite common in extended-range forecasts of 2-m temperature from ECMWF (Bertossa et al. 2021, hereafter B21). While this previous study identified the presence of univariate bimodality in forecasts, it did not analyze its cause. The goal of this paper is to identify flow configurations that lead to

bimodal forecasts, as a step toward understanding the processes and conditions under which bimodality can arise. We do so by introducing a novel method to identify such configurations, classifying them as “bimodal events.” This methodology clearly demonstrates that bimodal events can be linked to a diverse set of well-documented spatially and temporally coherent weather events.

The impact of particular atmospheric phenomena on the behavior of ensemble spread is an important topic of study, since the error growth of an ensemble is known to vary with the state of the atmosphere (Leutbecher and Palmer 2008; Palmer 2000). A specific type of distribution that may arise from ensemble spread is a bimodal distribution. Bimodal distributions occur when some ensemble members in a probabilistic forecast spread in such a way that they form a distinct group separate from the rest of the ensemble. Depending on how a continuous probability density function (PDF) is determined from the underlying discrete ensemble (a process sometimes referred to as “dressing”; Bröcker and Smith 2008), this can result in a corresponding PDF that contains two modes

Corresponding author: Cameron Bertossa, bertossa@wisc.edu

<footnote>ª Department of Atmospheric and Oceanic Sciences, University of Wisconsin–Madison, Madison, Wisconsin</footnote>

<footnote>ᵇ Department of Earth and Atmospheric Sciences, Cornell University, Ithaca, New York</footnote>

<footnote>ᶜ LMD-IPSL, Ecole Polytechnique, Institute Polytechnique de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Paris, France</footnote>
(Vannitsem et al. 2018). The method used by B21 to identify univariate bimodality was shown to be robust to these choices in dressing methods.

Identifying the presence of bimodality can be important in practice because it changes the interpretation of the forecast relative to a Gaussian assumption. This can strongly affect the forecast probability of crossing “critical thresholds” for decision makers managing risks, improvements to the forecast skill through postprocessing methodology, and a better understanding of the potential attraction to different atmospheric states (B21). However, detecting bimodality requires dressing methods that are computationally expensive. Understanding the origins of events that cause bimodality can contribute an indication as to when and where higher-order dressing methods should be used.

In this study, the ensemble forecasts of 2-m temperature identified as bimodal following B21 are used to identify associated coherent temporal or spatial flow patterns. The following questions guide our analysis: 1) Can we identify bimodal events or is bimodality incoherent in space and time? If we can identify bimodal events, what length and time scales are they typically associated with? 2) Can we connect bimodal events with particular atmospheric phenomena? 3) Is there evidence that different processes may be leading to bimodal events in different seasons or geographical regions?

To answer these questions, individual bimodal forecasts are grouped together based on which ensemble members reside in each mode. This allows us to identify coherent patterns of bimodality across forecasts without explicit prior consideration of the spatial and temporal scales. Once these coherent forecast groups are found, we quantify their spatial and temporal properties to answer question 1. Then, since each grouping is defined by two coherent modes, one is able to (i) examine the different geophysical characteristics associated within each mode and (ii) understand what processes may have led to their development, both of which contribute to answering question 2. Finally, examining these properties as a function of season and location naturally addresses question 3. We anticipate that relevant processes leading to bimodality will differ from region to region and from season to season, and so our goal here is simply to demonstrate that bimodal events exist and to lay the basis for further investigations, not to exhaustively characterize all such events.

This manuscript is organized as follows: In section 2 the cluster analysis used to identify bimodal events is explained. In section 3 this methodology is applied to three regions of interest, the properties exhibited by each region are examined, and physical explanations are hypothesized. Finally, section 4 concludes.

## 2. Data and methodology

As in B21, the forecasts considered here are from the European Centre for Medium-Range Weather Forecasts (ECMWF) Atmospheric Model Ensemble extended forecast (ENS extended) ensemble. This dataset includes 46-day forecasts of 2-m temperature, for a 50-member ensemble, initialized every Monday and Thursday (Haiden et al. 2019) from 3 December 2015 to

28 January 2021. The extended-range forecasts are produced by continuing the medium-range ENS forecast, but at lower resolution (horizontal resolution of around 36 km). ENS is an ensemble of 51 forecasts with a horizontal resolution of around 18 km. It comprises one control forecast plus 50 forecasts each with slightly altered initial conditions and slightly altered model physics. Extended-range forecasts are coupled to the Wave model (ECWAM) and the Dynamic Ocean model (NEMO) run with resolution of around 28 km (Persson and Grazzini 2007).

The methodology of B21 is used to identify whether a univariate forecast (hereinafter, “forecast” refers to an ensemble as a whole, not a single member, predicting a univariate quantity at a single grid point and lead time) is bimodal or not. Once a forecast has been identified as bimodal, the methodology provides a best guess of which ensemble members lie in each mode. This is done by first dressing the ensemble with a kernel density estimate (KDE), creating a continuous PDF of the field of interest (this study focuses on 2-m temperature) from the discrete ensemble. Then, depending on where ensemble members lie relative to the fit KDE’s relative minimum, they are assigned to be in the cold mode or the warm mode of the forecast (see Fig. 1a for an example). These mode assignments are then used to define a binary membership vector $\mathbf{x} = (x_1, x_2, ..., x_{N_E})$, where $N_E$ is the number of members in the ensemble; specifically, if member $i$ of an ensemble is in the cold mode of the ensemble, $x_i = 1$; if it is in the warm mode $x_i = 0$ (see Fig. 1b for an example). This is done individually for every bimodal forecast originating from the same initialization date (i.e., across grid points and lead times).

A clustering routine (Sasirekha and Baby 2013) can then be used to find similarities between membership vectors and form groups of bimodal forecasts for each initialization date. The specific clustering routine, its parameters used for this study, and a more exhaustive description of the methodology can be found in the appendix material.

A limitation of this approach, however, proves to be its computational cost; the clustering routine scales as $O[n^2 \log(n)]$ where $n$ is the number of binary membership vectors, i.e., the number of grid points and lead times in the domain under consideration. We thus limited the scope of the analysis in several ways. First, the analysis is performed on a few select regions rather than the entire globe (Fig. 1c exhibits an example of one of these regions). Regions are chosen based on high occurrence rates and large average separation between modes as determined by B21. Regions are also chosen to sample varied geographical configurations and are characteristically affected by very different atmospheric phenomena. Second, the analysis is only performed for the second and third week of forecast lead times. These weeks exhibit the highest rates of bimodality occurrence in ECMWF forecasts, and it is over this period that we expect the bimodality to differ most in character from that of the climatological distribution. Finally, the analysis only focuses on the largest bimodal events by restricting attention to only the largest five clusters identified in each region, for each forecast initialization date.
<figure><img src="image_1.png" /><figcaption>FIG. 1. Methodology flowchart. Each panel includes a short synopsis and a section label indicating where the process is outlined. (a) Tick labels represent ensemble members. Colors indicate which mode the ensemble member is assigned to based on its position relative to the local minimum of a fit kernel density estimate. (b) The binary standardization procedure from the forecasts in (a) is depicted. The distance metric as defined in section 2 is computed for the two forecasts. (c) The evolution of a cluster and its fit ellipse are depicted. The left panel is the evolution of a cluster's fit ellipse area. The entirety of lead times the clustering routine considers is provided (weeks 2 and 3 of forecast lead times). The green triangle represents when a bimodal event would start; the red triangle is when that event would end. Blue × tick marks represent the two particular forecast lead times depicted for the cluster in the right panel. For the right panel, green ticks represent all forecasts that are bimodal at that given lead time. Blue tick marks represent those bimodal forecasts that belong to the same cluster. An example of the ellipse fitting for these two time steps is presented. The black arrow indicates the movement of the ellipse center from one time step to the next and can be used as a means to calculate the discrete speed (listed at the bottom) and direction of propagation.</figcaption></figure>
While these choices introduce selection bias and are ultimately arbitrary, this approach is sufficient to achieve the goals laid out in the introduction above. We are able to demonstrate below that coherent bimodal events do occur in this forecast system, and are able to infer some details of the character of the largest and most coherent bimodal events in these regions and at these lead times. Further refinement of the methodology will likely be required should a more comprehensive study be deemed worthwhile.

### a. Cluster definitions

In this section, several definitions to help characterize each cluster are presented. These terms are referenced throughout the study.

It will prove useful to characterize how consistently the ensemble can be divided into two modes for a given cluster. We consider the number of times a given ensemble member lies in the cold mode of forecasts in a given cluster relative to the total number of forecasts within the cluster. If there are $N_C$ forecasts (univariate ensembles) in a cluster and $\mathbf{x}^l$ is the membership vector of the forecast j, the fraction of forecasts $f_i$ for which member i is in the cold mode is given by

$$f_i = \frac{1}{N_C} \sum_j^{N_C} \mathbf{x}_i^j \quad (1)$$

Each member is considered to be coherent ($c_i$) if it lies in one mode or the other in at least 95% of the forecasts that comprise the cluster. That is,

$$c_i = \begin{cases} 0 & \text{if } 0.05 < f_i < 0.95 \\ 1 & \text{otherwise.} \end{cases} \quad (2)$$

We then define the coherency of a cluster to be the number of members in the ensemble that are coherent across the cluster:

$$C = \sum_{i=1}^{N_E} c_i. \quad (3)$$

To provide an example of this measure of coherency, a synthetic example cluster is plotted in Fig. 2. This is an idealized case where the “cluster” is simply two points in space over a portion of the forecast lead times. The forecast set (hereinafter, we use the term “forecast set” to broadly refer to a group of forecasts which are related through a common grid point or lead time) in Fig. 2a is bimodal from all leads 15 through 30 and the forecast set in Fig. 2b is bimodal from leads 20 through 35. In practice, clusters are found to be made up of many different locations (e.g., Fig. 1c). However, this figure offers a more intuitive example of coherency. The dashed red line in both forecast sets represents the same member, which has especially low mode consistency in the cluster. For a portion of the bimodal lead times, this particular member is in the cooler mode of forecast set a but in the warmer mode of forecast set b. Furthermore, in both forecast sets, the ensemble member crosses over the relative minimum and resides in the opposite mode for some time. Both of these factors

<figure><img src="image_2.png" /><figcaption>FIG. 2. A synthetic example to depict coherency. Each panel represents a different forecast ensemble’s evolution as a function of lead time. Forecasts are made up of 50-member ensembles. Both forecast sets (set a and set b) belong to a single cluster that contains only these two spatial points but for the entirety of each forecast sets’ bimodal lead times. Gray shading indicates when an ensemble is bimodal. Thick black lines indicate the fit KDE’s relative minima. Dark blue lines indicate ensemble members that have high coherency ($f_i < 0.05$ or $f_i > 0.95$) whereas redder lines indicate low coherency. Correspondingly colored lines in forecast set a and forecast set b represent the same ensemble member. Dashed and bolded lines are selected members discussed in the text.</figcaption></figure>

contribute to its poor mode consistency within the cluster as a whole. In contrast, the dark blue member with the yellow outline represents a member with high mode consistency. It resides in the warmer mode of both forecast sets for the entirety of their bimodal periods. This particular cluster has a coherency of $C = 48$ members, meaning 48 out of the 50 members have values of $f_i > 0.95$ or $< 0.05$.

We define the representative modes of each cluster as the ensemble grouping that is most representative of the cluster. This is the membership vector (refer back to Fig. 1) that occurs most frequently within the forecasts belonging to this cluster; put another way, it is the membership vector with the highest probability of being drawn when sampling vectors randomly from the cluster. The representative modes are realized by individual ensembles more frequently for clusters that have very high coherency.

Note that we make a distinction between the occupancy and size of a cluster. The occupancy refers to the total number of univariate bimodal forecasts at all grid points and lead times within the cluster (in the particular case for Fig. 2, the occupancy is 32 since there are two grid points, each with 16 bimodal lead times). The size of a cluster refers to its
<figure><img src="image_3.png" /><figcaption>FIG. 3. (a) From Bertossa et al. (2021), the occurrence of bimodal forecasts in ECMWF 50-member 2-m temperature forecasts for November–December–January–February–March (NDJFM). (b) As in (a), but for March–June–July–August–September (MJJAS). (c) The distribution of occupancy for the largest (in occupancy) five clusters generated from each ECMWF forecast from the cold season of each region (NDJFM for the Northern Hemisphere and MJJAS for the Southern Hemisphere). Colors align with the boxed regions in (a) and (b). (e) Histogram of the number of members in each ensemble with at least 95% coherency ($f_i > 0.95$ or $<0.05$) for the clusters generated from the cold season of each region. (d) As in (c), but for the warm season of each region (MJJAS for the Northern Hemisphere and NDJFM for the Southern Hemisphere). (f) As in (e), but for the warm season of each region.</figcaption></figure>

spatial extent at a single lead time (or an average over a set of lead times). The size of the cluster in Fig. 2 may be very small, since it only consists of two forecast locations.

To further characterize the size and shape of a cluster, we fit an ellipse to the points within the cluster at every lead time. Explicitly, this is done using a covariance error ellipse (which takes advantage of the eigenvectors to find the directions of maximum variance; Amen 2012), where 90% of the cluster points at a given lead time are contained within the confidence ellipse. The area of this ellipse is then used as the representative area of the cluster at that lead time.

Several other properties of the cluster may be characterized with such an ellipse. We use the ellipse center as a representative location of the cluster. From the movement of the center from lead time to lead time, we can estimate the velocity of the cluster. An example of this process for a pair of lead times is depicted in the right panel of Fig. 1c. In this case, the ellipse, and thus the cluster, propagates eastward with speed of roughly 18 m s⁻¹.

b. Cluster behavior

The three particular regions that are studied in this manuscript are outlined in Figs. 3a and 3b. Figures 3c–f depict the occupancy and coherency of the largest five clusters found in each forecast for the three regions’ cold and warm season.

Each region's clusters are generally larger (in terms of occupancy) in the cold season than in the warm season. The South American region generally has the smallest clusters, while the Southern Ocean region has the largest clusters. Cluster coherency is generally worse for clusters with larger occupancy than clusters with lower occupancy (not shown); likely tied to this property, all three regions exhibit greater coherency in their respective warm seasons as compared to their cold seasons. The South American region exhibits the greatest coherency, while the Southern Ocean region exhibits the least. That being said, all three regions have a large population of clusters with a coherency of at least 25 members.

The spatial and temporal properties of the clusters that form may be assessed by plotting how the average cluster changes as a function of lead time. Since, however, the leads
times over which once cluster exists have no relation to those
for which another may be present, a reference lead time must
be chosen for each cluster to meaningfully compare their tem-
poral evolution. The reference lead time for each cluster is
chosen to be the lead time at which the cluster occupancy is
largest (lag 0 in Fig. 4). Then, for each cluster, the lead times
are aligned relative to the reference lead time and the occu-
pancy is normalized by the occupancy at the reference lead
time. The result of this analysis for each region is plotted in
Figs. 4a and 4b.

One property that emerges for all three regions is that clus-
ters, on average, monotonically increase in occupancy, reach
a peak, and then monotonically decrease as a function of lead
time. This temporal dependence is in no way implied by the
definition of the clusters; rather, this property emerges natu-
rally from the analysis.

Similar to occupancy, the coherency of the average cluster
grows steadily, reaches a peak, and then decreases sur-
rounding the lead time in which the occupancy is greatest;
though, there is the exception of a small dip in coherency at
a lag of 0 (Figs. 4c,d). Generally, this stands in contrast to
the relationship exhibited in Fig. 3, in which coherency de-
creases as cluster occupancy increases. This suggests that
these clusters are developing around meteorological events
that are evolving in systematically distinct ways in the differ-
ent ensemble members.

Clusters are frequently found to be concentrated together
in space, with one example being that from Fig. 1c. This
implies that the ellipse fitting methodology is a reasonable
approximation for the cluster area. As with occupancy and
coherency, the size of the typical cluster appears to steadily
grow, reach a peak, and then decrease (Figs. 4e,f). The
North Atlantic and South American regions have similar-
sized events on average, especially in the summer, whereas
the Southern Ocean's events are typically smaller. There may,
however, be some amount of bias due to the effect of latitude
on gridbox size. All three regions have clusters that are on the
order of 10⁵⁻¹⁰⁶ km² for several days, with the average cluster
being larger (in size) during the cold season versus the warm
season.

The combination of Figs. 4a–f supports the existence of co-
herent spatial and temporal structures that are resolved with
this clustering methodology. To reiterate, this clustering ap-
proach by design does not know about the spatial and tempo-
ral relationships between the forecast grid points that have
been identified as bimodal by the univariate methodology of
B21. While the coherency and occupancy of the clusters are
likely to be sensitive to the exact parameters of the clustering
algorithm used (e.g., one could likely find clusters with larger
occupancy if lower coherency was allowed), the tendency to
identify spatially and temporally coherent structures strongly
suggests a meteorological origin. In other words, this cluster-
ing approach has identified what could reasonably be termed
bimodal events, answering in the affirmative the first key sci-
ence question identified in the introduction. Clusters not only
appear to capture the maturing and deformation of an event,
but also have significant coherency, where well over half of

the ensemble members are consistently in one mode but not
the other for the majority of forecasts.

c. *Defining a bimodal event*

Supported by indications of cluster spatial and temporal
coherence (Fig. 4), we choose to refine each cluster into a bi-
modal event. We define bimodal events by only the lead times
in which the cluster occupancy continuously exceeds 10% of
its maximum value. If a cluster exceeds this threshold, dips
below the threshold and then re-exceeds the threshold, the
longest continuous stretch of lead times is used. This allows
one to define an event's persistence, or length, based on the
longest continuous length of lead times in which this threshold
is satisfied. The use of this threshold not only allows us to
cleanly define one bimodal event per cluster, but also reduces
noise in the spatial properties of the bimodal event that may
be present when only a few forecast points make up the clus-
ter (also making the ellipse representation more effective).
An example of the size evolution for a single cluster is pre-
sented in Fig. 1c. In this particular case, there is only one lead
time in which the cluster is defined but the bimodal event is
not (around 22 January).

For each region, the bimodal events derived from the five
largest clusters include about 10% of the total bimodal fore-
cast occurrences during weeks 2 and 3 of forecast lead times
found in Bertossa et al. (2021). It is clear, at least with how
the clustering process has been defined for this study, that
many of the total bimodal forecasts are not being used. Part
of this (likely a large portion) is a function of the clustering al-
gorithm and the suboptimal parameters chosen for it; much of
the signal of these bimodal forecasts has been thrown out
along the way by the choices made in this study. See, for ex-
ample, the larger extent of bimodality that appears to be spa-
tially coherent, yet is not included within the cluster in Fig. 1c.
Though unnecessary for answering the guiding questions pre-
sented here (as will be shown), the clustering procedure
should be refined in the future to take full advantage of the bi-
modal forecast set and reduce biases imparted by only exam-
ining the largest bimodal events.

**3. Regional analysis**

Thus far this study has shown the existence of bimodal
events, spatially and temporally coherent outbreaks of bi-
modality. Furthermore, tools (such as ellipse fitting) and
metrics (coherency, size, persistence, velocity) with which
these events can be studied and compared have been pre-
sented. Whether these tools that have been developed in
this study can help explain what meteorological phenomena
lead to the emergence of bimodal events is now explored. In
this section, a combination of ellipse properties, spatial oc-
currence maps and case studies are used to give plausible
explanations for phenomena resulting in bimodal events for
each of the three regions of interest. Each presented case
study is chosen such that it represents a majority of the cases
analyzed for each region. A summary section follows with
similarities between regions.
<figure><img src="image_4.png" /><figcaption>FIG. 4. (a) The mean cluster occupancy as a function of lead time for each region's cold season. Cluster sizes are normalized based on the maximum cluster size at a single lead time. Lags indicate the lead time preceding or following when the maximum cluster occupancy occurred. (b) As in (a), but for the warm season. (c) As in (a), but the y axis is the cluster coherency at each lead time. (d) As in (a), but the y axis is the cluster area according to a fit covariance error ellipse. (e) As in (a), and (c), but for each region's warm season. Some of the clusters included in this study may reach a maximum size very near the time bounds of our cluster analysis (beginning of week 2, end of week 3). Thus, in order to avoid inaccurate representations of evolution near these boundaries, the forecasts whose occupancy maxima occurred within the first 2.5 days or last 2.5 days of our analysis were not included in this figure (approximately 30% of the cases were removed).</figcaption></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 5. (a) PDF (normalized histogram) of bimodal event length for South America's cold season (MJJAS). Event length is defined by the longest continuous stretch of cluster lead times that exceed occupancies greater than 10% of the maximum occupancy. (c) Peak bimodal event area based on the lead time in which the area of the fit ellipse is the greatest. (e) Bimodal event velocity based on the movement of the ellipse center for each bimodal event. Note that a single event's velocity is determined based on the median direction of motion. Height of bar indicates the relative proportion of the direction of propagation for all bimodal events. Bars point toward direction of propagation. For a given direction, colors show the relative proportion of events that propagate at that speed. The mean speed of all events is listed in the legend. (b),(d),(f) As in (a), (c), and (e), but for the warm season (NDJFM).</figcaption></figure>

### a. South American region

The analysis of the South American region begins with a short discussion of the distributions of event length, event area, and event velocity based on fit ellipses for the region's bimodal events. These properties are presented for the cold season in Figs. 5a, 5c, and 5e, and for the warm season in Figs. 5b, 5d, and 5f. In general, events persist for roughly 1–2 days with slightly longer events occurring on average during the cold season as compared to the warm season. In agreement with Fig. 4, events are on the order of $10^5–10^6$ km$^2$. However, there is a noticeable size difference between the seasons, with the average cold season event area being twice as large as the average warm season event. Events in this region propagate at a speed of approximately $10–15$ m s$^{-1}$ on average, with slightly faster speeds being more common during the cold season. The bimodal events' direction of motion is typically eastward, with less frequent occurrences of northwest/southeast propagation, which becomes slightly more prevalent during the cold season.

Figure 6 depicts where bimodal events typically occur within the region, based on the number of forecast days for each season in which a bimodal event is present. Figure 6 suggests that bimodal events in the warm season commonly occur off the eastern coast of Brazil and Argentina and propagate eastward across the Atlantic. During the cold season, an additional set of bimodal events is identified around the eastern edge of the Andes, reaching a peak rate of nearly 40 forecast days per season. By analyzing individual tracks it can also be seen that the events near the Andes appear to be relatively stationary, with slight movement from south to north. This is distinct compared to the continuous eastward propagation associated with the Atlantic maximum. The alignment along the Andes may explain the increased northwest-southeast propagation seen in Fig. 5e.

To explore the Andes maximum more thoroughly, a specific bimodal event is analyzed in Fig. 7. This forecast set is initialized on 10 August 2017, during the region's cold season, and is valid for 19 August 2017. The properties of each mode
Bimodal Event Density

<figure><img src="image_11.png" /></figure>

FIG. 6. (a) South American bimodal event density (black and gray shaded), bimodal event start (blue contour) and end (red contour) locations, and two randomly sampled individual tracks (yellow and orange) depicting the direction of movement (arrows) for forecasts initialized during the warm season. The number of individual tracks plotted is heavily thinned from the total number of bimodal events (hundreds) for clarity. Event start (end) contours begin at 4 event start (end) locations per extended season within a 5° × 5° box, with an additional contour every 2 counts. Shading indicates the number of event days per season (e.g., the number of forecast days identified within a bimodal event; see the appendix). (b) As in (a), but for cold season initializations.

of this event are analyzed by taking the mean of the ensemble members that have been assigned to one mode or the other based on the cluster's representative modes. This particular event has a coherency of 27 members.

<figure><img src="image_12.png" /></figure>

FIG. 7. An ECMWF forecast set initialized on 10 Aug 2017 and valid for 19 Aug 2017. (a) Black contours represent the mean 850-hPa geopotential height contour for the ensemble members assigned to mode 1. Arrows indicate the mean 10-m winds for the ensemble members assigned to mode 1. Shaded contours indicate the mean 2-m temperature of the ensemble members assigned to mode 1. (b) Black contours represent the mean 850-hPa geopotential height contour for the ensemble members assigned to mode 2. Blue contours represent the difference in the mean 850-hPa geopotential height between mode 2 and mode 1. Arrows indicate the difference in 10-m winds between mode 2 and mode 1. Shaded contours indicate the difference in 2-m temperature between mode 2 and mode 1. The black ellipse indicates how the clustering algorithm has fit the event. Black contour intervals are 30 m. Blue contour intervals are 25 m.

Notable features in the geopotential height field in mode 1 (Fig. 7a) include the high pressure systems to the west and east of the continent which are separated by a trough with an axis roughly around 50°W. The particular location of the western high in Fig. 7a causes it to interact with the topography of the Andes mountains, and the eastern edge of the high wraps around the Andes. The large majority of the accompanying winds to these isobars are geostrophic with some ageostrophic flow toward the northern portion of the event. This brings cold air from the south as can be seen in the temperature field.

Figure 7b displays a succinct summary of mode 2, as well as the differences relative to mode 1. Mode 2 lacks almost any wavelike behavior in the geopotential heights and isobars remain relatively zonal. The west-east high-low-high departure in mode 1 is quite evident in the blue contours in Fig. 7b. The
<figure><img src="image_13.png" /><figcaption>FIG. 8. As in Fig. 5, but for the Southern Ocean region.</figcaption></figure>

lack of the wrapping high around the Andes leads to much warmer central South America temperatures in mode 2 because cold air advection from high latitudes is not occurring.

The geopotential height and wind structure exhibited in mode 1 of Fig. 7 is characteristic of a cold air incursion. This phenomenon is described in detail in Garreaud (2000). These events lead to cold air damming on the eastern edge of the Andes, qualitatively agreeing with the pattern evident in Fig. 7a. See Garreaud (2000, their Fig. 4b) for a comparison. Cold surges similar to the cold air intrusion event exhibited in Fig. 7 have been linked to Rossby wave breaking, a nonlinear atmospheric phenomenon (Sprenger et al. 2013). Not only do many of the bimodal event cases examined for this study mimic the two modes of Fig. 7, but these two states also represent the two most common modes of variability according to principal component analysis conducted in Compagnucci and Salles (1997). This may support the presence of regime-like behavior in this region to which the ensembles may be attracted (Michelangeli et al. 1995).

## b. Southern Ocean region

Figure 8 depicts the event persistence, area and velocity for the Southern Ocean region's cold and warm season. Typically, bimodal events in this region persist for roughly 2 days. The

typical event length is slightly longer for the cold season as compared to the warm season. Additionally, the average event area in the cold season is over double that of the warm season. Events in the Southern Ocean region propagate at a mean speed of 5–10 m s⁻¹, with faster events occurring on average during the cold season. The direction of propagation is mainly from west to east in the cold season, with a smaller fraction of propagation to the southeast. During the warm season, the direction of propagation is much more evenly split, with some southward and southwestward motion as well.

The event density and individual tracks for the Southern Ocean's cold and warm season are plotted in Fig. 9. In the cold season, Southern Ocean events mainly propagate eastward, north of the Antarctic Peninsula. However, during the warm season, bimodal events are more concentrated near the Ronne Ice Shelf and events do not appear to travel as far.

We present the temperature of each mode for all bimodal events within the Southern Ocean region (Fig. 10). Figure 10 indicates that most of the time, bimodal events' warm mode is very near the freezing point, irrespective of the time of year. In contrast, the cooler of the two modes varies greatly depending on when the forecast is initialized, with periodic behavior that aligns with seasonality. The lack of seasonal variability in the warmer mode that is very clearly evident in the
Bimodal Event Density

<figure><img src="image_14.png" /><figcaption>FIG. 9. As in Fig. 6, but for the Southern Ocean region. (a) Event start (end) contours begin at 2 event start (end) locations per extended season, with an additional contour every 2 counts. (b) Event start (end) contours begin at 3 event start (end) locations per extended season, with an additional contour every 3 counts.</figcaption></figure>

cooler mode points toward a constant boundary forcing that exists in the warmer mode but not the cooler. As of 2016, ECMWF included a dynamic-thermodynamic sea ice model as a part of ensemble forecasts, allowing for changing boundary conditions as the forecast evolves (Buizza et al. 2017). The fact that the warmer mode is almost always near the freezing point (close to typical sea surface temperatures in the polar oceans) while the cooler mode is well below freezing, indicates that the warmer mode is representative of a state with open ocean while the cooler mode corresponds to an ice-covered state. In the former, the relatively warm ocean can give away heat to the overlying atmosphere, helping regulate the temperature irrespective of the atmospheric state (also resulting in little spread in the warm PDF of Fig. 10b). In contrast, if this communication with the ocean is cut off (e.g., sea ice being present), 2-m temperature values have greater variability based on day-to-day atmospheric conditions (hence the greater spread in the cold PDF of Fig. 10b).

Figure 11 depicts a forecast set initialized during the Southern Ocean cold season, valid for 3 July 2019. Mode 1 has a weak low to the west of the Antarctic Peninsula. Mode 2, in

contrast, has a ridge axis in roughly this same location with a low further downstream. Extreme temperature departures between the modes occur both to the west and to the east of the peninsula, aligning with the departures in low-level winds. Accompanied with these temperature departures are differences in the sea ice fields, with reduced sea ice along the ice edge to the west of the peninsula for mode 2, and increased sea ice to the east. Notably, much of the bimodality in this region (whether that be a part of this particular bimodal event or not) develops near the ice edge.

Using this as a guide, the bimodality within this region could be developing in several different ways: differing advective influences, where flow from ocean-ice versus ice-ocean will have very different influences on the local 2-m temperature; or, by changes in the extent of the ice itself (e.g., melting or advection of sea ice).

The latter case is discussed first. In the context of Fig. 11, mode 2 has warm air advection from the north which may push (or melt) ice southward, exposing the relatively warm ocean, resulting in warmer 2-m temperatures in the western portion of the map. Similarly, cold air advection from the
<figure><img src="image_15.png" /><figcaption>FIG. 10. The most common temperatures of each mode for every bimodal event in the Southern Ocean region during MJJAS and NDJFM. Each blue (red) tick represents the most common temperature value for all ensemble members that belong to the cold (warm) mode in the bimodal forecasts of a particular bimodal event. A gray line has been marked for 0°C. (b) A PDF representation of the points presented in (a).</figcaption></figure>

south may advect (or form) sea ice further north, removing the temperature flux from ocean to atmosphere in these northerly locations, resulting in colder temperatures in the eastern portion of the map. Should this surface change occur in one mode but not the other, bimodality would certainly form. In the case of Fig. 11, bimodality is notably present between mode 1 and mode 2's ice edge contour.

That being said, the majority of bimodal points do not occur only where the sea ice fields between the two modes depart; in some cases, bimodality occurs quite far away from the ice edge. In these instances, the difference in the direction of low-level winds over the ice-ocean interface has the dominant role in developing the bimodal 2-m temperatures. Evidently, there is enough thermal inertia in the airmasses being advected from the ocean to the ice or vice versa to lead to the phenomenon exhibited in Fig. 10 (i.e., a near-freezing mode). The time scale of these bimodal events approximately aligns with previous studies of airmass “persistence” (i.e., how long an airmass can maintain its characteristics when displaced to a different location; Kotas et al. 2013).

The interaction of low-level winds with ice (whether that be sea ice or the ice sheet) appears to be the main cause of bimodality in this region. In the case of advection over sea ice, nonlinearity can be introduced if atmospheric waves deepen due to the sea ice induced temperature gradient. Several previous studies have connected Rossby wave breaking events to skewness in the climatological distribution of 850-hPa temperatures surrounding the southern hemispheric storm track (approximately 38° and 60°S) (Tamarin-Brodsky et al. 2019; Garfinkel and Harnik 2017). These two regions align with local maxima

in bimodality as apparent in Fig. 3b, with a minimum in between that may be representative of the center of the storm track axis; thus, similar processes may lead to skew versus bimodality. Alternatively, the development of preferential low-level wind modes may be attributed to blocking type regimes caused by the Antarctic Peninsula (Schwerdtfeger 1975; Massom et al. 2008; Hosking et al. 2013; Elvidge et al. 2016; Laffin et al. 2021).

### c. North Atlantic region

Figure 12 depicts histograms of event persistence, area, and velocity for the North Atlantic region's cold and warm season. The typical event persists for approximately 2 days, with slightly longer events occurring on average during the cold season. As seen in the previous two regions, the typical event area approximately doubles in size for forecasts initialized during the cold season versus during the warm season. Events propagate at an average speed of 5–10 m s⁻¹, with slightly faster events occurring on average during the cold season. Roughly two-thirds of the time, events have an eastward direction of propagation in the cold season, with the remainder of the motion being predominantly southeastward. However, the direction of propagation is more evenly split between eastward and southeastward motion for warm season forecasts.

The event density and individual tracks for the warm and cold season for the North Atlantic region are plotted in Fig. 13. During the summer season, events are concentrated between 70°–35°W and 35°–45°N, with a maxima of over 25 event days per season. During the North Atlantic cold season, bimodal events appear to occur over a greater range of locations, with
<figure><img src="image_16.png" /></figure>

FIG. 11. The two modes of a bimodal event from a forecast set initialized on 17 Jun 2019, and valid for 3 Jul 2019. (a) The sea level pressure (black contour) contoured every 4 hPa, 10-m winds (arrow), and 50% sea ice concentration line (blue contour) for mode 1. Shaded contours indicate the difference in the sea ice concentration between mode 2 and mode 1. Green ticks represent all bimodal points at this lead time. (b) As in (a), but for mode 2. Shaded contours now indicate the difference in 2-m temperature between mode 2 and mode 1. Blue ticks represent only those bimodal points that belong to the bimodal event at this lead time.

relatively even frequencies over much of the region depicted; however, there are notable maxima located near the Gulf of Mexico and the Gulf Stream. Individual tracks mainly progress west to east, with some northwest to southeast propagation especially prevalent in the cold season's Gulf Stream maximum. The persistent Atlantic maximum in Fig. 13 may indicate that sea surface temperature gradients near the Gulf Stream play a dominant role in the evolution of bimodal events in this region. However, given that these events propagate at a speed of 5–10 m s⁻¹, this suggests that events are not fixed to the SST gradients, but rather, atmospheric processes play a significant role in their evolution as well.

Figure 14 depicts the evolution of a forecast set initialized during the Northern Hemisphere cold season, valid for 10 and 11 February 2016. The bimodality in this forecast is associated with the displacement of a persistent large-scale high pressure system that is present in the two modes. Mode 1 is characterized by the development of a smaller low pressure system to the east of the high, whereas this is absent in mode 2. The high in mode 1 is displaced west (possibly due to the

formation of the low), relative to the high in mode 2. This leads to warm anomalies as the high in mode 2 has northerly flow over the equator-to-pole SST gradient.

This persistent subtropical high was a major factor in most of the bimodal events analyzed in this region. If the high is displaced westward enough, it leads to anomalously warm conditions in the southern part of the United States as air is advected northward from the Gulf of Mexico.

Persistent sea surface temperature gradients also appear to play a significant role, allowing for departures between the low-level winds in each mode to induce large differences in the 2-m temperature. The displacement or deformation of a subtropical high seems to be the main cause for the circulation differences in the bimodal events examined here. Such a system directly affects the low-level wind's interaction with the SST induced gradient. What causes this displacement or deformation is unclear and likely variable. The North Atlantic region has been identified as an area for frequent Rossby wave breaking (Homeyer and Bowman 2013) as well as bomb cyclones, rapidly intensifying low pressure systems (Sanders
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 12. As in Fig. 5, but for the North Atlantic cold (NDJFM) and warm (MJJAS) season.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

and Gyakum 1980). These systems quickly develop due to their interaction with Gulf Stream. As such, they may introduce the nonlinearity required for multiple modes to form in an ensemble.

That being said, as previously mentioned, there are several distinct event locations within this region, such as southern Canada and the Gulf of Mexico. As such, the Gulf Stream is clearly not the only forcing mechanism leading to bimodality in this region. Alternate meteorological phenomena which may lead to the development of bimodality within this region include polar lows that develop near Hudson Bay (Gachon et al. 2003).

### d. Summary discussion

While the character of the bimodal events identified in each region is distinct, surface features that favor the advection of distinctly different air masses come out as a frequent ingredient in the cases examined here: a north-south topographic barrier, or a strong temperature gradient in the underlying surface such as the sea ice edge or a western boundary current. That being said, the typical propagation speeds found in each region suggests that the bimodality in 2-m temperature is not fixed relative to the underlying features.

In the Southern Ocean, the distinct development of synoptic-scale waves in each mode resembles similar processes in

the North Atlantic in which strong sea surface temperature gradients play an important role. While the South American region is distinct in the sense that its bimodal events are not directly connected to a temperature gradient, the boundary conditions (via topography) still emerge as an important ingredient. Furthermore, in some Southern Ocean cases, topographic blocking due to the Antarctic Peninsula may mirror the role that the Andes have in the South American region.

In all three regions, bimodal events can be connected to impactful atmospheric phenomena, whose distinct character would likely be underpredicted should a Gaussian distribution be assumed. In the case of the North Atlantic, the displacement of the Bermuda high is an important topic of study, since it has drastic effects on temperature and precipitation (Li et al. 2011; Zhu and Liang 2013; Diem 2013). In the South American region, cold surges during the winter can lead to flash freezing of crops, whereas summertime episodes can influence rainfall (Garreaud 2000). Finally, bimodal events in the Southern Ocean region may be connected to changes in the sea ice around Antarctica or the ice sheet itself.

### 4. Conclusions

This paper introduces a novel methodology to identify and characterize “bimodal events,” that is, meteorological conditions
Bimodal Event Density

<figure><img src="image_23.png" /><figcaption>FIG. 13. As in Fig. 6, but for the North Atlantic region. (a) Event start (end) contours begin at 2 event start (end) locations per extended season, with an additional contour every 1 counts. (b) Event start (end) contours begin at 1 event start (end) locations per extended season, with an additional contour every 2 counts.</figcaption></figure>

that give rise to spatially and temporally coherent bimodality in subseasonal-to-seasonal 2-m temperature forecasts from ECMWF. While there are still pathways for refinement of the approach presented here, the premise has proven useful in demonstrating that the univariate bimodality identified by Bertossa et al. (2021) is in fact connected to coherent atmospheric phenomena, a primary goal of this study. Bimodal forecasts are found to have spatial and temporal coherence that can be identified solely by which ensemble members are present in each mode of the univariate bimodal forecasts. While being able to characterize events based on similarities in forecast behavior makes sense in theory, it is remarkable how well the idea behaves in the context of bimodality. The authors are unaware of such an approach being used previously to characterize weather events, despite the practicality of wanting to understand how specific forecast behavior can be linked to weather phenomena.

The computational expense of the methodology restricted our considerations to three specific regions: South America, the Southern Ocean, and the North Atlantic. On average, the bimodal events identified here last 1–3 days and cover a spatial area on the order of $10^5–10^6$ km². However, there are

distinctions between particular events, between different regions, and between seasons. For instance, within the South American region, bimodal events that occur near the Andes have tracks that are relatively stationary with some south-north propagation, while in the western Atlantic, bimodal events tend to propagate steadily eastward at 5–10 m s⁻¹.

We further examined the meteorological development of the two modes in specific case studies in each region (Figs. 7, 11, and 14). The combination of general characteristics and case studies of bimodal events indicates several weather phenomena that contribute to the development of coherent bimodality in forecast ensembles. In the South American region, the most common bimodal event appears to be due to one mode developing a blocking pattern which results in cold air incursions along the eastern flank of the Andes, while the other mode has flow that remains relatively zonal. In the Southern Ocean, low-level winds' interaction with sea ice appears to be the main cause of bimodality in forecasts. Finally, in the North Atlantic, bimodal events are associated with the displacement or deformation of a persistent subtropical high.

While the specific processes leading to bimodal events in each region appear to be distinct, each are related to interactions between large-scale dynamical processes with more local geographic features at the surface including topography, sea surface temperature gradients, and the margins of the cryosphere.

We have suggested hypotheses for processes that are nonlinear in nature, and thus can possibly lead to the development of two modes. However, it is difficult to confirm these without a much more extended analysis. The fact that differences between distinct subsets of forecast ensembles can grow while members of each subset remain similar suggests the presence of distinct feedback mechanisms within each mode, but the process (or processes) responsible may vary from region to region or even from event to event. To answer these questions, more systematic composite studies or the use of conceptual models are needed. Our results suggest several processes that such conceptual models might need to include.

Undoubtedly, given the computational limitations of the present approach and the exploratory nature of our methodology, we have not exhaustively revealed all processes that lead to the development of bimodality in ensemble forecasts. With further improvements or a more geographically targeted approach, this methodology may be adapted to characterize bimodal events in a more complete manner, uncovering similar insights. Several potential improvements include a more computationally efficient clustering algorithm that would allow for analyzing larger domains, considering different geophysical variables, extending the clustering algorithm across multiple initialization dates to understand if bimodal events continue to persist closer to the validation date, and revising the normalization factor and distance thresholds to capture a larger fraction of the univariate bimodality.

Nonetheless, our approach has proven useful for identifying coherent sets of ensemble members that represent different modes, something that has not yet been done. This has allowed us to identify several common ingredients that are associated with the development of large-scale bimodal events
<figure><img src="image_24.png" /><figcaption>FIG. 14. The two modes of a bimodal event from a forecast set initialized on 25 Jan 2016, and valid for 10 and 11 Feb 2016. Mode 1 contains 5 members, and mode 2 contains 45 members. (left) The sea level pressure (black contour) contoured every 2 hPa, 10-m winds (arrow), and sea surface temperature (shaded) for mode 1. (right) As in the left panels, but for mode 2. Shaded contours now indicate the difference in 2-m temperature between mode 2 and mode 1. The lead time that each row of panels is valid for is listed at the top right of each panel. In all panels, blue ticks represent the bimodal event locations at that lead time.</figcaption></figure>

(i.e., large-scale atmospheric processes interacting with local boundary conditions). Accurately representing these processes in models and capturing their initial conditions is thus vital to achieving skillful forecasts. Moreover, we have been able to connect bimodal events to previously identified specific weather phenomena. There is no indication in the literature that these particular phenomena should lead to bimodal or, arguably, “regime-like” behavior in ensemble forecasts;

thus, this work presents a new perspective on these phenomena. A deeper understanding of the dynamics and the prediction of these meteorological phenomena may be gained by further study of their effects on forecast distributions.

*Acknowledgments.* We thank Cornell University for funding and the use of materials and facilities, as well as the European Centre for Medium-Range Weather Forecasts (ECMWF) for
providing the dataset used. The authors also acknowledge the data center ESPRI/IPSL for their help in storing and accessing the data.

*Data availability statement.* The ECMWF ENS extended forecasts are accessible through the ECMWF website (https://www.ecmwf.int/en/forecasts). ERA-Interim data can be found at http://www.ecmwf.int/research/era.

# APPENDIX

## Additional Methodology

### a. Clustering routine

Clusters of similar forecasts based on ensemble membership vectors are formed using a hierarchical clustering routine (Sasirekha and Baby 2013), specifically Python's scipy.cluster. hierarchy.fclusterdata function (scipy v1.4.1) which performs hierarchical clustering using the single linkage algorithm with a custom distance formula.

The custom distance formula used to determine the similarity between two forecasts' membership vector **x** and **y** is defined as

$$d(\mathbf{x}, \mathbf{y}) = \frac{1}{\sqrt{C_x C_y + W_x W_y}} \sum_{i=1}^{N_E} (x_i - y_i)^2, \quad (A1)$$

where $N_E$ is the number of members in the ensembles, $C_x = \sum_i x_i$ is the number of members in the cold mode of ensemble **x** and $W_x = N_E - C_x$ is the number of members in the warm mode of ensemble **x**. The terms $C_y$ and $W_y$ are defined similarly. The expression is symmetric, in the sense that $d(\mathbf{x}, \mathbf{y}) = d(\mathbf{y}, \mathbf{x})$. The normalization factor is largest when the cold modes of both forecasts are either much smaller or much larger than the warm modes ($C_x \ll W_x$ and $C_y \ll W_y$, or $C_x \gg W_x$ and $C_y \gg W_y$). Conversely it is smallest when the cold modes of each forecast are very different sizes ($C_x \gg C_y$ and $W_y \gg W_x$, or $C_x \ll C_y$ and $W_y \ll W_x$). It has been included to reduce the impact of a single member being in a different mode when one of the two modes is much larger than the other in both forecasts. Other normalization factors are also possible and several other expressions were also tested; the form (A1) is found to give the best results (generate the largest clusters with lower intracluster distances). An example of this process for two synthetic 5-member forecasts is depicted in Fig. 1b.

What determines a cluster is an iterative routine involving a cluster center, whose distance is minimized from surrounding points. The number of clusters can be predetermined either by a priori information or via a “distance threshold,” defined as $d_c$. The threshold $d_c$ specifies the maximum intercluster distance allowed, where smaller distances would imply more clusters; in this case, the clustering algorithm decides how many clusters will be present. As the threshold approaches 0, each point would be its own cluster. The value for $d_c$ used in this study is 1/37. This value was chosen through trial and error based on maximizing cluster coherency and occupancy.

This is a similar method to which the normalization factor for Eq. (A1) was tested.

Note that clustering is performed on each set of initialized forecasts (where forecasts are initialized twice a week) independently. That is, the clustering routine only searches for similar binary membership vectors for forecasts with identical initialization dates (which are valid across various grid points and lead times). Since only the largest five clusters are considered from each initialization date, and there are approximately 5 years' worth of forecast data (with two initialization per week), we have a sample size of approximately $2 \times 52 \times 5 \times 5 = 2600$ clusters for each region. Each of these clusters may include bimodal forecasts from any grid point within the region's domain, as well as any lead time of week two and three.

### b. Calculation of bimodal event density

Bimodal event density (Figs. 6, 9, and 13) is calculated based on where ellipse centers lie for each lead time that a bimodal event exists. Note that different forecast initialization dates may have a bimodal event at the same location for the same validation date, since validation dates can overlap for different forecast initialization. That is to say, for example, 30 event days for an extended season that is 150 days long does not necessarily mean that 20% of the days are predicted to be bimodal events, but rather, of all the forecast lead times analyzed within a given season, 30 days' worth, on average, have a bimodal event occurring. Since forecast sets are initialized twice per week, 14 days per forecast are considered (weeks 2 and 3), and there are approximately 21 weeks, this equates to around 588 forecast days per extended season.

# REFERENCES

Amen, S. K., 2012: Linear estimation for data with error ellipses. Ph.D. thesis, The University of Texas at Austin, 48 pp., https://repositories.lib.utexas.edu/handle/2152/ETD-UT-2012-05-5270.

Bertossa, C., P. Hitchcock, A. DeGaetano, and R. Plougonven, 2021: Bimodality in ensemble forecasts of 2 m temperature: Identification. *Wea. Climate Dyn.*, **2**, 1209–1224, https://doi.org/10.5194/wcd-2-1209-2021.

Bröcker, J., and L. A. Smith, 2008: From ensemble forecasts to predictive distribution functions. *Tellus*, **60A**, 663–678, https://doi.org/10.1111/j.1600-0870.2008.00333.x.

Buizza, R., J.-R. Bidlot, M. Janousek, S. Keeley, K. Mogensen, and D. Richardson, 2017: New IFS cycle brings sea-ice coupling and higher ocean resolution. *ECMWF Newsletter*, No. 150, ECMWF, Reading, United Kingdom, 14–17, https://www.ecmwf.int/en/newsletter/150/meteorology/new-ifs-cycle-brings-sea-ice-coupling-and-higher-ocean-resolution#top.

Compagnucci, R. H., and M. A. Salles, 1997: Surface pressure patterns during the year over southern South America. *Int. J. Climatol.*, **17**, 635–653, https://doi.org/10.1002/(SICI)1097-0088 (199705)17:6<635::AID-JOC81>3.0.CO;2-B.

Diem, J. E., 2013: Influences of the Bermuda high and atmospheric moistening on changes in summer rainfall in the Atlanta, Georgia region, USA. *Int. J. Climatol.*, **33**, 160–172, https://doi.org/10.1002/joc.3421.
Elvidge, A. D., I. A. Renfrew, J. C. King, A. Orr, and T. A. Lachlan-Cope, 2016: Foehn warming distributions in non-linear and linear flow regimes: A focus on the Antarctic Peninsula. *Quart. J. Roy. Meteor. Soc.*, **142**, 618–631, https://doi.org/10.1002/qj.2489.

Gachon, P., R. Laprise, P. Zwack, and F. J. Saucier, 2003: The effects of interactions between surface forcings in the development of a model-simulated polar low in Hudson Bay. *Tellus*, **55A**, 61–87, https://doi.org/10.3402/tellusa.v55i1.12079.

Garfinkel, C. I., and N. Harnik, 2017: The non-Gaussianity and spatial asymmetry of temperature extremes relative to the storm track: The role of horizontal advection. *J. Climate*, **30**, 445–464, https://doi.org/10.1175/JCLI-D-15-0806.1.

Garreaud, R., 2000: Cold air incursions over subtropical South America: Mean structure and dynamics. *Mon. Wea. Rev.*, **128**, 2544–2559, https://doi.org/10.1175/1520-0493(2000)128<2544:CAIOSS>2.0.CO;2.

Haiden, T., M. Janousek, F. Vitart, L. Ferranti, and F. Prates, 2019: Evaluation of ECMWF forecasts, including the 2019 upgrade. ECMWF Tech. Memo. 853, 56 pp., https://www.ecmwf.int/node/19277.

Homeyer, C. R., and K. P. Bowman, 2013: Rossby wave breaking and transport between the tropics and extratropics above the subtropical jet. *J. Atmos. Sci.*, **70**, 607–626, https://doi.org/10.1175/JAS-D-12-0198.1.

Hosking, J. S., A. Orr, G. J. Marshall, J. Turner, and T. Phillips, 2013: The influence of the Amundsen–Bellingshausen Seas low on the climate of West Antarctica and its representation in coupled climate model simulations. *J. Climate*, **26**, 6633–6648, https://doi.org/10.1175/JCLI-D-12-00813.1.

Kotas, P., R. Twardosz, and Z. Nieckarz, 2013: Variability of air mass occurrence in southern Poland (1951–2010). *Theor. Appl. Climatol.*, **114**, 615–623, https://doi.org/10.1007/s00704-013-0861-9.

Laffin, M. K., C. S. Zender, S. Singh, J. M. Van Wessem, C. J. P. P. Smeets, and C. Reijmer, 2021: Climatology and evolution of the Antarctic Peninsula föhn wind-induced melt regime from 1979–2018. *J. Geophys. Res. Atmos.*, **126**, e2020JD033682, https://doi.org/10.1029/2020JD033682.

Leutbecher, M., and T. N. Palmer, 2008: Ensemble forecasting. *J. Comput. Phys.*, **227**, 3515–3539, https://doi.org/10.1016/j.jcp.2007.02.014.

Li, W., L. Li, R. Fu, Y. Deng, and H. Wang, 2011: Changes to the North Atlantic subtropical high and its role in the

intensification of summer rainfall variability in the southeastern United States. *J. Climate*, **24**, 1499–1506, https://doi.org/10.1175/2010JCLI3829.1.

Massom, R. A., S. E. Stammerjohn, W. Lefebvre, S. A. Harangozo, N. Adams, T. A. Scambos, M. J. Pook, and C. Fowler, 2008: West Antarctic Peninsula Sea ice in 2005: Extreme ice compaction and ice edge retreat due to strong anomaly with respect to climate. *J. Geophys. Res.*, **113**, C02S20, https://doi.org/10.1029/2007JC004239.

Michelangeli, P.-A., R. Vautard, and B. Legras, 1995: Weather regimes: Recurrence and quasi stationarity. *J. Atmos. Sci.*, **52**, 1237–1256, https://doi.org/10.1175/1520-0469(1995)052<1237:WARRAQS>2.0.CO;2.

Palmer, T. N., 2000: Predicting uncertainty in forecasts of weather and climate. *Rep. Prog. Phys.*, **63**, 71–116, https://doi.org/10.1088/0034-4885/63/2/201.

Persson, A., and F. Grazzini, 2007: User guide to ECMWF forecast products. ECMWF Meteor. Bull. M3.2, ECMWF, 162 pp.

Sanders, F., and J. R. Gyakum, 1980: Synoptic-dynamic climatology of the “bomb.” *Mon. Wea. Rev.*, **108**, 1589–1606, https://doi.org/10.1175/1520-0493(1980)108<1589:SDCOT>2.0.CO;2.

Sasirekha, K., and P. Baby, 2013: Agglomerative hierarchical clustering algorithm—A review. *Int. J. Sci. Res. Publ.*, **3**, 1–3.

Schwerdtfeger, W., 1975: The effect of the Antarctic Peninsula on the temperature regime of the Weddell Sea. *Mon. Wea. Rev.*, **103**, 45–51, https://doi.org/10.1175/1520-0493(1975)103<0045:TEOTAP>2.0.CO;2.

Sprenger, M., O. Martius, and J. Arnold, 2013: Cold surge episodes over southeastern Brazil—A potential vorticity perspective. *Int. J. Climatol.*, **33**, 2758–2767, https://doi.org/10.1002/joc.3618.

Tamarin-Brodsky, T., K. Hodges, B. J. Hoskins, and T. G. Shepherd, 2019: A dynamical perspective on atmospheric temperature variability and its response to climate change. *J. Climate*, **32**, 1707–1724, https://doi.org/10.1175/JCLI-D-18-0462.1.

Toth, Z., O. Talagrand, G. Candille, and Y. Zhu, 2003: Probability and ensemble forecasts. *Forecast Verification: A Practitioner's Guide in Atmospheric Science*, I. T. Jolliffe and D. B. Stephenson, Eds., John Wiley and Sons, 137–164.

Vannitsem, S., D. S. Wilks, and J. Messner, 2018: Statistical Post-processing of Ensemble Forecasts. Elsevier, 320 pp.

Zhu, J., and X.-Z. Liang, 2013: Impacts of the Bermuda high on regional climate and ozone over the United States. *J. Climate*, **26**, 1018–1032, https://doi.org/10.1175/JCLI-D-12-00168.1.