ERA5‐Based Significant Tornado Environments in Canada Between 1980 and 2020

This study uses ERA5 close‐proximity soundings and associated convective parameters to characterize significant tornadic storm (F/EF2+) environments between 1980 and 2020 in parts of Canada. It is shown that ERA5 convective parameters are suitable to represent observed parameters, based on radiosonde comparisons. Results indicate that the eastern Canadian Prairies have nearly double the lifted condensation level with higher level of free convection compared to eastern Canada (southern Ontario/Quebec). Eastern Canada has more a humid boundary layer and free troposphere that can lead to warmer cold pools, favoring tornadogenesis. Central Canada (Manitoba) has the largest mixed‐layer (ML) convective available potential energy (CAPE) mainly due to a combination of regional differences in low level moisture and steeper mid‐level lapse rates in western Canada. Central continental U.S. and Canadian regions appear to have the highest (most negative) convective inhibition, leading to more explosive initiation. Mean bulk wind shear and storm relative helicity (SRH) increases from west to east, with eastern regions being significantly larger. The supercell composite and significant tornado parameters are generally less than U.S. magnitudes, particularly in western Canada, and would require recalibration for more practical use in Canada. Overall, western Canada significant tornadic storms are associated with more low‐level looping hodographs and dominated by thermodynamic influences compared to larger wind influences in eastern regions. This is likely due to more spring, late summer, and autumn events that typically have well‐developed synoptic systems (stronger wind shear) with overall less buoyant energy in eastern regions.


Introduction
On average in any given year, Canada can experience 60 (based on actual reports) to 230 (based on modeling) tornadoes (Cheng et al., 2013;Newark, 1984;Sills et al., 2012), resulting in multi-millions of dollars in damage and loss of life (e.g., Sills et al., 2012).The Barrie 1985 (Etkin et al., 2001) and Edmonton 1987 (Bullas & Wallace, 1988) tornadoes are examples of the most costly and significant tornado events in Canadian history.
Despite their impacts in Canada, there have been limited studies to better understand tornadic storm environments, primarily due to a lack of data.Most published work focused on Alberta due to their prolific hailstorms, and southern Ontario with the highest population density.Alberta research focused on hailstorms and general conceptual models of severe thunderstorm environments (e.g., Chisholm & Rennick, 1972;Strong, 1979Strong, , 1986;;Smith & Yau, 1993a, 1993b), however, Dupilka and Reuter (2006a, 2006b, 2011) summarized various convective parameters for tornadic storms in central Alberta, using one radiosonde site.They suggested the degree of lowlevel veering, along with storm relative helicity (SRH), were associated with stronger tornadic storms; thresholds were 0-3 km SRH > 150 m 2 s 2 and 900-500 mb shear exceeding 3 m s 1 km 1 for F2-F4 tornadoes, which are smaller than typical U.S. environments (e.g., Thompson et al., 2003).
In Ontario, a relation between lake breezes and tornado climatology was inferred with limited data (King, 1997;King et al., 1996); Ontario lake breezes can initiate convection, but may also play a role in tornadogenesis via lake breeze-storm interactions (e.g., King et al., 2003).Sills and King (2000) looked at landspout tornadoes generated by coincident misovortices and growing convective towers along lake breeze boundaries, however, cases were limited and only applied to non-mesocyclone events.
Most tornadic storm mesoscale environments have been studied outside of Canada, using close-proximity model or reanalysis soundings.The reader is referred to Coffer et al. (2020) and Taszarek et al. (2020) for a review with only limited discussion here, focusing on tornadic storms only.Examples of early work include Brooks et al. (2003) and Thompson et al. (2003) who used model and reanalysis-derived convective parameters to characterize tornadic environments.Shortly after, newer model and reanalyzes datasets were used to further advance knowledge, make comparisons amongst U.S. regions and convective modes (e.g., Anderson-Frey et al., 2016, 2019;Coffer et al., 2019;Gensini et al., 2021;Grams et al., 2012;Thompson et al., 2007Thompson et al., , 2013)).All of these studies pointed to the importance of various convective available potential energy (CAPE) and wind shear layers, including low level SRH, with recent work suggesting entrainment CAPE may lead to more accurate updraft velocities (Peters et al., 2020).For example, a drier middle troposphere could promote more entrainment within the updraft and weaken it, while stronger storm-relative inflow may enhance it (Peters et al., 2020).Taszarek et al. (2020) (hereafter T2020) compared European and U.S. convective environments, including tornadoes, using ERA5 (Hersbach et al., 2020) and found that the U.S. has higher moisture, CAPE, convective inhibition (CIN), wind shear, and mid-tropospheric lapse rates, whereas Europe has higher 0-3-km CAPE and low-level lapse rates.They also noted that the supercell composite parameter (SCP) and significant tornado parameter (STP) were calibrated over the U.S. Coffer et al. (2020) also used ERA5 to show that 0-500 m SRH can be useful for tornado prediction in the U.S. and Europe, however, 100-200 m SRH layers were the most skillful in Europe.A recent ERA5 study comparing China and U.S. strong tornado (F/EF2 + ) events suggested that China has less favorable kinematic conditions than the U.S., which partially explains its less frequent strong tornado events (Zhang et al., 2023).The objectives of this article are to: (a) quantify important tornado environment convective parameters associated with significant tornado (F/EF2 + ) events using ERA5 in regions of Canada between 1980 and 2020, and (b) statistically distinguish Canadian regional significant tornado environments.Brief comparisons to U.S., Europe and China are also made where possible to do so.The broader goal is to provide knowledge of "typical" significant tornado environments, and their variations, to contribute to prediction capabilities.Various ERA5 vertical profile-derived convective parameters that are important to describe tornadic environments are estimated, similar to Anderson-Frey et al. (2016), Gensini andBrooks (2018), , andGrams et al. (2012).However, we will closely follow T2020 since this study uses a similar, but expanded, ERA5-derived convective parameter dataset that were processed with the same computational script.This article is organized as follows: Section 2 presents the data and methods, including the study area, various datasets used, statistical analysis methods, and limitations of the work.Section 3 contains the main results, focusing on a brief ERA5 convective parameter validation, then regional tornadic storm environments for various thermodynamic, kinematic and composite parameters, as well as composite thermodynamic and wind profiles.Section 4 concludes by summarizing the key findings and future research.

Tornado Data and Study Area
This study used F/EF2 + tornado events from the Canadian tornado database between 1980 and 2009 (Sills et al., 2012) (https://open.canada.ca/data/en/dataset/fd3355a7-ae34-4df7-b477-07306182db69)and 2010-2020 data made available by the Northern Tornadoes Project (NTP; Sills et al., 2020).As of early 2024, some NTP source data from Environment and Climate Change Canada between 2010 and 2016 are not yet public.Relevant to this study, data included: (a) date and time of tornado, (b) nearest town/city and province, (c) approximate latitude/longitude of witnessed tornado during its lifecycle, (d) F/EF rating.Canada began using the EF-scale in April 2013 (Sills et al., 2014) and was implemented in the database after this date.Although difficult to quantify, errors associated with the F/EF scale transition are not expected to be significant for F/EF2 + events; <2% increases (shifts) in F2 + from F to EF scales (Edwards et al., 2021).
Most of the significant tornadoes (F/EF2 + ) occurred between Alberta and Quebec, thus, the focus is on those 5 provinces.Ontario was divided into northern (west of 85°W) and southern (east of 85°W) regions due to: (a) its vast size longitudinally, where storm environments may be different, and (b) the tornado "gap" north and northeast of Lake Superior (Figure 1 and Cheng et al., 2013).Abbreviated region names are used hereafter (Figure 1, Table 1).Subdividing other provinces was not done due to limited tornado samples that can affect statistical comparisons.Grouping tornado events by provincial regions was done to simplify the analysis, but still allows for reasonable environment assessments based on the authors' experience.Future work will examine cluster analysis or machine learning to group Canadian tornado events.There were 211F/EF2 + tornadoes between 1980 and 2020 with 35 days of multiple (two or more) significant tornadoes; several cases did not have sufficient information (e.g., date and/or time) and were discarded.For the 35 days with multi-tornado reports, we minimized environmental statistical oversampling by applying a "Goldilocks zone" (GZ) for optimal proximity soundings to be within 1-2 hr and 40-80 km of the actual tornado based on Potvin et al. (2010).Therefore, only the highest F/EF-rated occurrence was kept in the analysis for any multi-tornado day, unless other tornadoes were >80 km or more than 2 hr from the highest F/EF-rated event.If two or more tornadoes had the equivalent highest-rated rank on any day, and occurred within the GZ of one another, then only one of those tornadoes was randomly selected to "represent" the environment.After applying these criteria, a total of 166 F/EF2 + cases with unique ERA5 profiles were used in the analysis (45 removed) (Table 1; Figure 1).84% were F/EF2, while 12% were F/EF3, with only 4% F/ EF4 cases (Table 1).The only recorded F/EF5 in Canadian history occurred on 22 June 2007 in Manitoba.Between 1980 and 2020, N_ON (MB) did not have any recorded F/EF2 + cases prior to 1994 (1984), while SK and AB first cases were 1982 and 1981, respectively; causes for these may be due to population density changes over time that can affect reports or problematic reports that were filtered out of the national database.
Seasonal and diurnal distributions of strong tornado events are important to highlight since this can help explain regional environment differences.Eastern Canadian (S_ON/QC) events are distributed over a longer seasonal period compared to western regions (Figure 2).In addition, S_ON and QC had 63% and 48% of their events occur in the April to May and August to  September periods, respectively (Figure 2).However, 87% of AB, 84% of SK and 63% of MB events occurred between June and July (Figure 2).If S_ON/QC are grouped together and compared to western regions (AB/SK/ MB combined), the diurnal peak in tornado occurrence occurs near 1645 solar time in the west compared to 1515 solar time in the east (Figure 3); the N_ON peak lies between eastern and western Canada (not shown).Solar times were used to mimic the diurnal heating period.The p-value (0.067) from a Kolmogorov-Smirnov test is suggestive of storms occurring 1.5 hr later in the day in western Canada compared to their eastern counterparts.

ERA5 Convective Parameters
Similar ERA5 data (Hersbach et al., 2020) and derived convective parameters were used as in T2020, so an abbreviated discussion is provided.The rationale for using ERA5 is that: (a) ERA5 represents North American and European convective environments reasonably well (Taszarek et al., 2021;hereafter T2021), and (b) direct comparisons between other ERA5 studies can be made.Important ERA5 aspects include (Hersbach et al., 2020): (a) 0.25°(31 km) horizontal resolution, (b) 137 terrain-following hybrid-sigma levels including 28 levels in the lowest 2 km for improved boundary layer representation, and (c) hourly temporal resolution.
ERA5 vertical profiles were collected for each of the 166 F/EF2 + tornadoes as close as spatially and temporally as possible to actual events (time = 0; T0), including ± 4 hr around T0; 9 profiles per tornado event (T-4, T-3, …, T0, …, T+3, T+4).Manual profile inspection was performed for each to ensure temporal errors were minimized; for example, if the ERA5 timing of convective precipitation (storm) was not accurate, this would potentially affect any profile and lead to its contamination by simulated precipitation (King & Kennedy, 2019).The "representative" pre-convective profile was selected based on boundary layer evolution (e.g., combination of maximum near-surface temperature timing and/or any obvious signs of precipitation or cold pool contamination) between T-4 to T+4 of each tornado event (Brooks et al., 1994).This led to 90% of the representative profiles to be within ±2 hr of tornado occurrence (most were T-2, T-1 or T0), with no T-4 or T+4 profiles used.Multiple profiles (i.e., grid points) in space were not considered since the selection of representative profiles in time would potentially alleviate some of the ERA5 spatial errors, and averaging several grid-point profiles to represent the environment (instead of one single profile) could overly smooth the profile.The drawback of not explicitly accounting for ERA5 spatial errors is that one assumes ERA5 storms are at least within a  reasonable distance (e.g., 80 km) of the actual storm.It is unknown if this criterion was met for every event.
This study used thundeR v.1.1 (Taszarek et al., 2023) to generate the convective and tornado parameters, similar to T2020 (Table S1 in Supporting Information S1).Only some parameters showed Canadian regional differences and are highlighted here.
Although T2021 compared ERA5 and observed soundings over North America, only limited Canadian detail was available.Thus, we conduct a brief ERA5 convective parameter comparison between operational soundings and ERA5 in Section 3.1; parameter definitions are shown in Table S1 in Supporting Information S1.The comparison is not meant to be exhaustive, but nonetheless, is important to ensure convective parameter confidence.There are only four Canadian sounding sites within the domain of interest (see Figure 1).Moosonee, ON was not included due to its proximity to James Bay that could pose local environment variations, and since no strong tornadoes occurred in this area.
Comparisons were made between convective parameters derived from the four operational sounding sites and nearest ERA5 grid point soundings, similar to T2021.Sounding data between 1990 and 2020 from the University of Wyoming were used (http://weather.uwyo.edu/upperair/).Stringent quality-control was applied to the operational soundings to remove incomplete, unphysical or error-prone soundings following the same methodology as T2021.Only time periods between 1800 and 0000 UTC May 1 and September 30 were investigated, to coincide with the convective season and typical peak heating hours.As a result, 16,005 rawinsonde measurements were compared with ERA5 profiles over the 31 year period.Several selected convective parameters and statistical error metrics were calculated for the operational and ERA5 soundings (Figures S1-S3 and Tables S2-S3 in Supporting Information S1).The mixed-layer LCL appear to have discrete levels at higher altitudes (Figure S1b in Supporting Information S1) that is caused by the thundeR script using a 50 m vertical resolution at higher LCLs to limit computational time.Locally estimated scatterplot smoothing (LOESS) lines were included to illustrate local deviations between observed and ERA5 parameters, similar to T2021 (Jacoby, 2000).Where there are sparse data points, the LOESS may be less accurate in depicting biases.All comparisons included operational soundings that had 0-500 m mixed-layer CAPE (ML_CAPE) > 0 J kg 1 .Emphasis was placed on ML parcels, as opposed to most-unstable (MU) parcels, since T2021 found better correlations for ML parcels; results here were similar in most cases (not shown).

Statistical Analysis for Regional Comparisons
To explore regional (geographic) differences of the individual convective parameters, univariate tests and visualizations were used.Most convective parameters were normally distributed, although some had heavy-tails and outliers.Given this, heteroscedastic one-way Analysis of Variance (ANOVA) for trimmed means were performed as the global test for the convective parameters (Wilcox, 2016) and Yuen's trimmed means for the pairwise post-hoc tests (Yuen, 1974).Heteroscedastic one-way ANOVA for trimmed means is a robust extension of one-way ANOVA using trimmed data (extreme values removed) that accounts for unequal variances across groups (Wilcox, 2016).Similarly, Yuen's trimmed means is a robust alternative to the t-test where extreme values are trimmed and where equality of variances is violated (Yuen, 1974).Because the post-hoc Yuen's trimmed means tests were used for comparing multiple combinations of regions, the P-values for these were adjusted using the Holm-Bonferroni method for multiple comparisons (Holm, 1979).
Using the above methods, univariate tests were summarized with boxplots (Figures 4-9) in R (Patil, 2021).A post-hoc comparison between two regions is indicated by a horizontal bracket connecting the groups.Horizontal brackets are used to "connect" regions that are statistically different (significant at an adjusted P ≤ 0.05), while convective parameters that are not significant had no connecting brackets.
To examine the multivariate relationships, Multiple Discriminant Analysis (MDA) was performed on the regionally classified convective parameters (Friendly & Fox, 2021) in CRAN R. MDA finds the linear combination of variables (ERA5 convective parameters) that best separate three or more groups (i.e., the regional locations), generating canonical discriminant axes (Wilks, 2011;dos Santos et al., 2023).Because MDA is a linear method, the convective parameters were first standardized using the Ordered Quantile Normalization procedure (Peterson, 2021;Peterson & Cavanaugh, 2020).We tested the significance of the canonical axes using Wilk's Lambda (Friendly & Fox, 2021) and used the significant axes to construct a biplot.Scores for each site on each canonical axis are plotted along with structure correlations for the original variables.The resulting biplot (Figure 10) was used to visually assess whether the geographical regions were statistically different (P ≤ 0.05) along with the influence of ERA5 parameters.Pairwise post-hoc t-tests (adjusted using the Holm-Bonferroni method) were used to assess statistically different geographic regions for the significant axes.Further information about MDA and how it is being applied to this work can be found in Supporting Information.

Limitations
Limitations of this work include: (a) Limited sample sizes in each region.The robust statistical tests provide greater confidence in regional comparisons, and it is believed there are enough meteorologically representative cases in each region, although more would be valuable, (b) The use of provincial boundaries to group tornado events.Future research will examine cluster analysis or machine learning techniques, nonetheless, this study did show statistically significant regional differences, (c) All reanalyzes have biases, timing and spatial errors in storm initiation, intensity and placement, and mesoscale effects, such as storm outflows, that may not be captured (e.g., Allen & Karoly, 2014;Gensini et al., 2014;King & Kennedy, 2019;Taszarek et al., 2018).This study attempted to account for timing (and indirectly spatial) errors, however, more rigorous analysis would be required to address this issue.S1 in Supporting Information S1.Convective variables are calculated from the ERA5 proximity grid point.

Results
First, we begin with a brief comparison between convective parameters derived from ERA5 and observed soundings in Section 3.1.Sections 3.2 onward addresses our primary objectives.Only some of the convective parameters are shown for brevity (Figures 4-9) that indicated the most statistical differences amongst regions, with Table S4 in Supporting Information S1 containing statistics for all 38 parameters.However, some important parameters such as convective inhibition were also shown, that did not indicate statistically different regions.

ERA5 Versus Observed Convective Parameter Comparison
ERA5-derived convective parameters replicated observed parameters well in most cases, with similar errors as T2021.For example, mixed-layer mixing ratio and dew points had correlations between 0.88 and 0.96 with low mean absolute errors (MAEs) and small positive biases (Figure S1 and Table S2 in Supporting Information S1).
Mixed-layer LCL, mixed-layer LFC, and mixed-layer CAPE were also well correlated to observations (0.65-0.93) with some negative biases and similar magnitude MAEs compared to T2021 (Figure S1 and Table S2 in Supporting Information S1).ERA5 appears to underestimate larger CAPE, shown by the LOESS, that was also noted by T2021.
ERA5 temperature parameters had high correlations with observations (all >0.95) (Figure S1 and Table S2 in Supporting Information S1) with small negative biases and slightly larger MAEs in low levels (10 m-850 hPa),  then improved at 700 and 500 hPa (Table S2 in Supporting Information S1).All lapse rates had correlations >0.9 except for 0-1 km (0.73).Results are similar to T2021.
Wind speeds were well correlated with observations, improving with height from 0.74 at 10 m to 0.97 at 500 hPa, with small negative biases and MAEs (Figure S2 and Table S3 in Supporting Information S1).Wind shear correlations mimicked wind speeds, with increasing correlations from low levels to higher levels, with small biases and MAEs all near 2 m s 1 (Figure S2 and Table S3 in Supporting Information S1).Finally, ERA5 SRH between 0 and 1 km and 0-3 km both showed good correlation to observations (0.82-0.85) with relatively small positive biases and MAEs (Figure S2 and Table S3 in Supporting Information S1).All wind parameter errors were similar to T2021.
Pointed out by T2021, composite parameters often suffer from the largest error due to the multiplying errors with each term; the more terms in a composite parameter, the larger the error.For example, SCP and STP (both consist of more than two thermodynamic and kinematic terms) had the lowest correlation (0.59-0.68) (Figure S3 and Table S3 in Supporting Information S1).However, the biases and MAEs are small (Table S3 in Supporting Information S1), and smaller than those in T2021.Both ERA5 composite parameters underestimate soundingderived parameters with increasing magnitude, as shown by the LOESS (Figure S3 in Supporting Information S1).Overall, it is not unexpected for the ERA5 validation to reveal acceptable comparisons to observed soundings since re-analyses assimilate these observations.This is particularly true in our study since the vast majority of comparisons are valid for standard times (0000 UTC).

Thermodynamic Parameters
Despite being situated in the center of Canada, MB has similar mixed-layer mixing ratios as eastern Canada, with SK and AB having lowest magnitudes (Figure 4a).With AB/SK being in the lee of the Rocky Mountains and void of large inland lakes, it is not surprising they are driest (e.g., Oke et al. (1998)), however, MB must have considerable moisture advection contributions.Apart from advection, the Canadian Prairies are reliant on local surface moisture sources, such as annual field crops (e.g., Hanesiak et al., 2004;Brimelow et al., 2011;Hanesiak et al., 2011) that have been linked to its tornado climatology (Raddatz & Cummine, 2003).AB mixed-layer mixing ratios are similar to strong AB tornado events based on observed soundings (Dupilka & Reuter, 2011).Canadian mixing ratio magnitudes are within those reported by T2020 for the U.S., with AB more closely aligned with Europe (T2020) and western U.S. (for example, Zipser & Golden, 1979;Szoke et al., 1984).
CAPE is important for updraft intensity, vertical stretching and therefore tornadic supercells (e.g., Brooks et al., 1994;Fawbush & Miller, 1954;Thompson et al., 2003Thompson et al., , 2004)).MB has the highest mixed-layer CAPE with QC being lowest (Figure 4b).Since the regional patterns of CAPE and mixed-layer mixing ratio are different, this suggests that low level moisture is not the only variable that determines regional CAPE differences across Canada.Grams et al. (2012) showed that U.S. northern plains events have higher CAPE (median 2100 ± 900 J kg 1 ) compared to other U.S. regions.SK/MB/N_ON fall into northern plains and central Chinese (Zhang et al., 2023) environments.In contrast, S_ON and QC CAPE are aligned with U.S. midwest environments.AB CAPE is comparable to Dupilka and Reuter (2006a) as well as Colorado case studies (e.g., Murdzek et al., 2020).Rasmussen (2003) suggested that CAPE between 0 and 3 km may be important for enhancing low level updrafts and vertical stretching in U.S. tornadic supercells.Although regional comparisons are not statistically significant here (due to variability) (Figure 4c), S_ON and QC have the largest mixed-layer 0-3 km CAPE and greater number of larger magnitude cases (dots in Figure 4c).Eastern Canada magnitudes are within the upper end of Rasmussen's U.S. and T2020 U.S. and European thresholds.
Convective inhibition (CIN) is critical to assess convection initiation likelihood, but also for the build-up of energy to produce explosive storm development (e.g., Fawbush & Miller, 1954;Maddox, 1976).Larger CIN can delay convective initiation to evening hours and result in more explosive storm development in a more favorable wind shear environment (wind shear can increase toward evening hours).In addition, CIN can also support isolated storm-mode, which is favorable to tornadoes and large hail.However, if CIN is too strong, it can prevent any storm development.Thus, sufficient, but not too strong CIN can support stronger, isolated storms.Although not regionally significant, MB/N_ON have the most negative mixed-layer CIN (Figure 4d), and comparisons with T2020, Brooks et al. (1994), Davies (2004) and Zhang et al. (2023) suggests that central North American environments may have higher convective inhibition compared to western and eastern counterparts, but similar to southern and central China.
Low mixed layer LCL and mixed-layer LFC have been linked to tornadic storms by lowering the mid-level updraft and increasing vertical stretching (e.g., Edwards & Thompson, 2000;Rasmussen & Blanchard, 1998;Thompson et al., 2003).Both mixed-layer LCL and LFC are generally higher in western Canada (Figures 5a and  5b).Eastern regions have the lowest LCL and LFC primarily due to more abundant boundary layer moisture and smaller dew point spreads (discussed in Section 3.6).Eastern Canada is near typical U.S. and European heights (e.g., Brooks et al., 1994;Davies, 2004;T2020;Pilguj et al., 2022), northern Ontario similar to China (Zhang et al., 2023), while western Canada is similar to U.S. northern plains (e.g., Thompson et al., 2003;Davies, 2004).
Lapse rates are intrinsically linked to associated CAPE, and they all suggested similar trends (Figures 5d and 5e).
Western regions were statistically more negative than eastern regions with 1-1.5 K km 1 steeper lapse rates.This in combination with the low level moisture contributes to the bell-shaped mixed-layer CAPE pattern centered over Manitoba.Eastern Canada has similar lapse rates as U.S., European and southern Chinese counterparts (e.g., Johns & Doswell, 1992;Rasmussen & Blanchard, 1998;T2020;Zhang et al., 2023).

Kinematic Parameters
Western (eastern) Canada mean wind magnitudes are smaller than (as large as) many U.S. strong tornado environments (e.g., Kerr & Darkow, 1996;Markowski et al., 2003).For example, the 1-3 km mean wind (used as a surrogate for the low-level jet (LLJ) in T2020) was highest in QC and dropped off markedly westward (Figure 6a).LLJ winds in western Canada are on the low end of European events, while eastern Canada are comparable to upper-end and median events in Europe and U.S., respectively (compared to Markowski et al., 2003 andT2020).
Wind shear is critical for convection organization, storm rotation and longevity (e.g., Johns & Doswell, 1992;Markowski et al., 1998).Bulk wind shear generally increased from west to east in Canada (Figures 6b-6e).Bulk shear between 0 and 1 km in Ontario is twice as large as western Canada, with Quebec having the largest values (Figure 6b).Bulk shear between 0-3 km and 0-6 km show similar trends (Figures 6c and 6d), while effective bulk shear has less range between west and east (Figure 6e).Western (eastern) Canada wind shear is smaller than (as large as) typical U.S. and European F/EF2+ tornado environments (e.g., Thompson et al., 2003;Markowski et al., 2003;T2020;Pilguj et al., 2022).China and western Canada have similar wind shear magnitudes (Zhang et al., 2023).MDA biplot with regional means, 95% confidence intervals (circles), and structure correlations of the convective parameters with the canonical axes (Can1 along the x-axis and Can2 along the y-axis).
Storm relative helicity (SRH) has been used to distinguish tornadic versus non-tornadic storms, with higher magnitudes linked to stronger tornadoes (e.g., Davies-Jones, 1984;Davies-Jones et al., 1990;Markowski et al., 2003;Thompson et al., 2003).SRH (for right moving storms) generally increases from west to east in Canada with western regions statistically smaller than eastern regions for all SRH layers (Figure 7).The 0-100 m and 0-500 m layers were included since Coffer et al. (2019Coffer et al. ( , 2020) ) showed they were the most skillful SRH parameters in predicting strong tornadic supercells, however, their median magnitudes (>80 m 2 s 2 for 0-100 m; >200 m 2 s 2 for 0-500 m) were than Canadian cases.They also noted that these lowest layers have better forecast skill in the south and eastern U.S. compared to the northern plains.Alberta 0-1 km SRH (0-3 km SRH) are slightly larger (smaller) compared to Dupilka and Reuter (2011), however, this could be due to different datasets as well as sample sizes.Western Canada 0-1 km SRH are smaller than the U.S. and northern/central Chinese counterparts, however, from Manitoba eastward, 0-3 km SRH are similar to U.S., European and Chinese magnitudes (compared to Coffer et al., 2020;Thompson et al., 2003;T2020;Zhang et al., 2023).Bunkers et al. (2000) storm motions were used to reveal mean (normal) storm tracks and right-mover deviant storm motions, since there has been no formal climatology of tornadic storm motions in Canada.However, this technique has biases, particularly for significantly tornadic supercells (Bunkers, 2018).Mean wind (normal) storm motion azimuth (Bunkers_MW_A) shifts from SSW in Alberta to WSW once in eastern Canada (Figure 8a) and increase in mean wind magnitude (Bunkers_MW_M) from west to east (Figure 8b).Results here show a more southerly component in Alberta storms compared to Dupilka and Reuter (2011) (250°at 13 m s 1 ), but similar velocity.The largest right-mover deviant storms occur in Manitoba, followed by Saskatchewan/Alberta, with eastern regions being similar (Figure 8c).Based on the authors' experience, the Bunkers mean and right-mover deviation storm track directions are reasonable, although individual storms/ cases can be different.

Other Indices and Composites
Other indices that showed west-east differences included the 0-2 km moisture flux (g s 1 m 2 ) and cold pool strength (K); see Table S1 in Supporting Information S1 for definitions.The moisture flux between 0 and 2 km signifies horizontal moisture flows that are critical for supplying storm energy, while the cold pool strength has been linked to tornadic versus non-tornadic supercells (e.g., Markowski et al., 2002).Large increases in 0-2 km moisture flux from west to east are noted with highest levels in S_ON and QC (Figure 9a).Based on the authors' experience and our low level moisture and mean low level wind results, these regional differences are not surprising since eastern Canada is a more humid climate (e.g., Oke et al., 1998).Larger cold pool strengths can lead to tornadogenesis failure (Markowski et al., 2002), and results suggest that western regions have significant tornadoes despite having colder cold pools (predicted by the method used here).Cold pool strength was highest in SK and MB and are statistically different from all eastern regions (Figure 9b).It should be noted that ERA5 profiles only provide larger scale pre-storm conditions, and likely do not capture internal storm-scale processes responsible for actual cold pools.
The supercell composite parameter (SCP) and significant tornado parameter (STP) are commonly used to assess supercell and tornado likelihood; see Table S1 in Supporting Information S1 (e.g., Coffer et al., 2019;Gropp & Davenport, 2018;Thompson et al., 2004Thompson et al., , 2007)).Only SCP from Gropp and Davenport (2018) and STP from Coffer et al. (2019) will be discussed, since patterns are similar to older SCP and STP formulations.SCP peaks in MB and decreases westward and eastward, with smallest values in AB (Figure 9c).The western Prairies are similar to typical northern plains, European and south China values (e.g., Grams et al., 2012;T2020;Zhang et al., 2023).Manitoba eastward are smaller than U.S. tornadic supercells in Thompson et al. (2003) but similar to north/central China (Zhang et al., 2023) and general U.S. supercells in Thompson et al. (2004).STP increases markedly from western to eastern Canada (Figure 9d) with western provinces similar to European and south China magnitudes and within the bottom end of U.S. northern plains (Grams et al., 2012) and north China (Zhang et al., 2023).Eastern Canada is near lower-end values of Thompson et al. (2004) and the midwest U.S., but similar to north/central China (Grams et al., 2012;Zhang et al., 2023).Overall, the STP, and possibly SCP, would have to be recalibrated for Canada, particularly in western Canada.Given that western Canadian events have less wind shear but higher CAPE than eastern Canada, we suspect adjustments to wind shear terms in the SCP and STP formulations may improve Canadian utility.

Regional Discrimination of Combined Parameters
To put regional comparisons from Sections 3.2-3.4into perspective and determine overall regional significance, MDA was performed to provide statistical tests and a biplot for interpretation (Figure 10); see Section 2.3.Unlike the univariate methods presented previously, this approach includes possible variable interactions that further discriminate among geographic regions.Family-wise error rate is also controlled as the method extracts composite discriminant axes that can be tested for group separation.The first two canonical discriminant axes were most significant (Wilk's Λ = 0.05, p-value < 2.2e 16 and Λ = 0.32, p-value = 1.6e 04 , respectively) while the third and subsequent discriminant axes (not shown) were not significant.
The convective parameters differ significantly amongst all regions in Canada as evidenced by the nonoverlapping confidence intervals (95% Figure 10) and confirmed with post-hoc pairwise tests (not shown).
The first discriminant axis (Can1 = 78.7%canonical discrimination) indicates that the most significant differences in convective parameters follow a west to east gradient.As one runs along the x-axis (Can1), the regional differences are significant (non-overlapping circles), except AB versus SK and S_ON versus QC (Figure 10).
There is significant but comparatively lesser regional separation along the second axis (y-axis Can2 = 9.4% canonical discrimination).Along the y-axis, the Prairie provinces are significantly different from each other (nonoverlapping circles).Quebec is different from all Ontario, AB is different than all Ontario, and MB is different from eastern regions.
In the top five positive convective parameters along the first axis, 0-1 km bulk shear has the largest structure correlation (0.69) with Bunkers mean wind magnitude, 1-3 km mean wind, 0-100 m SRH, and 0-500 m SRH having correlations >0.68 (Figure 10); these are all wind-related parameters that increase in magnitude eastward.Mixed-layer LCL, right-mover deviation (R_Mover_Dev) and cold pool strength had the largest negative structure correlations ( 0.62, 0.60 and 0.57 respectively) on the first axis and are strongly associated with western Canada.On the second axis, the largest positive structure correlation was 2-4 km lapse rate (0.19) but the highest absolute (negative) correlation was mixed-layer CAPE ( 0.55) followed by mixed-layer mixing ratio; most significant parameters are thermodynamic-related.

Composite Skew-T's and Hodographs
Composite ERA5 skew-t and hodographs were generated for each region (Figures 11 and 12).The skew-t "averaging" scheme stems from Warren et al. (2021) who used a vertical profile of relative humidity and parcel buoyancy for averaging.
The SK/MB skew-t's reveal "loaded gun" soundings with larger mixed-layer CAPE compared to other regions (Figures 4 and 11).However, the boundary layer is deeper and drier (i.e., larger dew point spreads) in SK/MB compared to eastern regions, leading to higher LCLs and LFCs (Figures 4 and 11) and higher convective temperatures (>30°C) compared to eastern regions (26-27°C).Interestingly, mixed-layer parcel convective cloud depths (i.e., mixed layer equilibrium level minus LFC) are largest in eastern Canada and smallest in Alberta, while eastern regions have skinnier mixed-layer CAPE compared to Saskatchewan/Manitoba (Figure 11).Convective cloud depth is larger in eastern Canada mainly due to lower LFCs there.Lower LCLs and LFCs in S_ON/QC also contributes to them having larger 0-3 km mixed-layer CAPE than any other region.Eastern Canada also has smaller dew point spreads above the boundary layer compared to Canadian Prairies, that can lead to less dilution of convective updrafts (e.g., Houston & Niyogi, 2007) and warmer cold pools (Markowski et al., 2002).However, eastern regions have less steep lapse rates (Figure 4) that can counter reduced dilution (Houston & Niyogi, 2007).Chinese composite skew-t's are most similar to eastern Canada.
Composite hodographs (Figure 12) reveal that western regions (particularly Alberta) have significantly more 0-1 km curvature and less bulk shear compared to eastern regions.Western regions have more of a directional and speed shear mix compared to eastern regions that have large speed shear and SRH in low levels (see also Figures 6  and 7).The shape of the AB hodograph is similar to Dupilka and Reuter (2011) and may be related to mountainplains processes (Strong, 1986), dryline effects (Taylor et al., 2011), or a combination of the two, mixed with general synoptic set-up (Smith & Yau, 1993a, 1993b).It appears that MB and China have the most similar hodographs.Interestingly, eastern Canada experiences larger 6-12 km winds (see Figure 11 wind barbs) that can lead to differences in storm speed, supercell type, and updraft width that can assist with higher mass fluxes to promote more precipitation and affect downdraft characteristics (e.g., Rasmussen & Straka, 1998;Warren et al., 2017).However, storm type and morphology are beyond the scope of this study.

Summary and Conclusions
This article uses 1980-2020 ERA5 reanalysis proximity soundings and several convective parameters (Table S1 in Supporting Information S1) to (a) quantify commonly used environment parameters associated with significant tornado (F/EF2 + ) events in parts of Canada, and (b) distinguish these regional environments using statistical tests.
Comparisons to other international studies are made where possible, to suggest future research of examining why similarities exist around the globe.ERA5-derived convective parameters represent observed parameters well, based on rawinsonde comparisons.After filtering the database, 166 F/EF2 + events remained in the analysis (Figure 1 & Table 1).A brief summary of the most notable results follows.
Western Canada (Prairies) has nearly double the LCLs and higher LFCs in eastern Prairies compared to eastern Canada (southern Ontario/Quebec) that has higher boundary layer moisture and smaller dew point spreads.In fact, eastern regions have smaller dew point spreads throughout the troposphere, that can lead to less convective updraft dilution (e.g., Houston & Niyogi, 2007) and warmer cold pools (Markowski et al., 2002) that favor tornadogenesis.However, internal storm-scale processes are largely responsible for cold pool temperature, which is not factored in our analysis.Despite higher moisture in the east, central Canada (Manitoba) has the largest mixed-layer CAPE, due to relatively high low-level moisture and steeper mid-level lapse rates in western Canada.Larger (and fatter) CAPE and generally larger convective inhibition in eastern Prairies can lead to more delayed, explosive storm development compared to eastern regions.In fact, peak occurrence of strong tornadoes in western regions occurs 1.5 hr later than eastern regions.This timing difference could be due to (a) greater capping in western regions, and/or (b) western region tornadoes developing later in the storm's life cycle.Overall, Chinese composite skew-t's are most similar to eastern Canada, while Saskatchewan/Manitoba are similar to northern U.S. Plains.
Mean bulk wind shear and low level SRH increases from west to east, with eastern regions being significantly larger, and can lead to differences in storm speed, supercell type, and updraft width that can assist with higher mass fluxes to promote more precipitation and affect downdraft characteristics (e.g., Rasmussen & Straka, 1998;Warren et al., 2017).Western regions have more looping hodographs with a directional and speed shear mix, compared to eastern regions that have large speed shear and SRH in low levels.Much larger low level SRH in eastern regions is primarily driven by larger speed shears.Overall, Manitoba has the most similar hodographs to China, while eastern Canada are most similar to eastern U.S. environments.
The SCP and STP are generally less than U.S. magnitudes, particularly western Prairies, and would require recalibration to be beneficial.Taszarek et al. (2020) also showed limited utility of current SCP and STP formulation for European environments.Since western Canada tends to have larger CAPE, we suspect that adjustments to wind shear terms within the SCP and STP formulations may improve the utility of these parameters in western Canada.Overall, canonical correlation analysis and MDA suggests that all 6 regions (Figure 1) are statistically distinct from each other, with respect to the 38 convective parameters used in this study.Of particular significance, MDA showed distinct west-to-east differences, however, further research is needed to explain the physical reasons for this.Nonetheless, it is shown that Canadian Prairie strong tornadic storms are more dominated by thermodynamic and buoyant energy influences compared to larger wind shear influences in eastern regions.We suspect eastern Canada has stronger dynamics and winds at all levels partially due to more Spring and Autumn events that typically have well-developed synoptic systems and over all less CAPE Future analysis will compare summer versus spring/autumn environments, since southern Ontario has been experiencing more late-season events (Sills et al., 2022) and future cool-season tornadoes may be stronger (Woods et al., 2023).

Figure 3 .
Figure 3.Diurnal probability density distribution of strong tornadoes (F/ EF2+) between 1980 and 2020 for eastern Canada (S_ON and QC combined) and western Canada (AB, SK, MB combined).The horizontal axis is mean solar time.

Figure 4 .
Figure 4. Box-and-whisker plots of (a) mixed-layer mixing ratio, (b) mixed-layer CAPE, (c) 0-3 km CAPE, and (d) mixedlayer CIN for each Canadian region.The median is the solid horizontal line inside the box, the box edges represent the 25th and 75th percentiles, whiskers represent the 10th and 90th percentiles, and circles are jittered raw data observations.Horizontal bars connecting the regions denote statistically significant differences.Convective parameter definitions are in TableS1in Supporting Information S1.Convective variables are calculated from the ERA5 proximity grid point.

Figure 10 .
Figure10.Summary of the canonical discrimination of all ERA5-derived convective parameters described in this paper.MDA biplot with regional means, 95% confidence intervals (circles), and structure correlations of the convective parameters with the canonical axes (Can1 along the x-axis and Can2 along the y-axis).