Hailstones are a natural hazard that pose a significant threat to property and are responsible for significant economic losses each year in the United States. Detailed understanding of their characteristics is essential to mitigate their impact. Identifying the dynamic and physical factors contributing to hail formation and hailstone sizes is of great importance to weather and climate prediction and policymakers. In this study, we have analyzed the temporal and spatial variabilities of severe hail occurrences over the U.S. southern Great Plains (SGP) states from 2004 to 2016 using two hail datasets: hail reports from the Storm Prediction Center and the newly developed radar-retrieved maximum expected size of hail (MESH). It is found that severe and significant severe hail occurrences have a considerable year-to-year temporal variability in the SGP region. The interannual variabilities have a strong correspondence with sea surface temperature anomalies over the northern Gulf of Mexico and there is no outlier. The year 2016 is identified as an outlier for the correlations with both El Niño–Southern Oscillation (ENSO) and aerosol loading. The correlations with ENSO and aerosol loading are not statistically robust to inclusion of the outlier 2016. Statistical analysis without the outlier 2016 shows that 1) aerosols that may be mainly from northern Mexico have the largest correlation with hail interannual variability among the three factors and 2) meteorological covariation does not significantly contribute to the high correlation. These analyses warrant further investigations of aerosol impacts on hail occurrence.
Severe convective storms with hail, tornadoes, and strong winds are responsible for significant property damages and economic losses in the United States each year. The economic losses from severe convective storms were $125 billion during 2007–16, more than all other natural hazards (i.e., flooding, droughts, and tropical cyclones) combined (Baggett et al. 2018). Most annual hail damages happen during the spring and summer months in the central United States (Changnon 1999; Cecil and Blankenship 2012; Allen et al. 2015b; Tippett et al. 2015). Hailstone damages to property are a function of kinetic energy, which can increase rapidly as the diameter increases. The hailstorm that occurred on 24 May 2011 in the Dallas–Fort Worth metroplex produced severe hail in excess of 4 in. (10.16 cm) in diameter and an estimated economic loss of $876.8 million (Brown et al. 2015).
Understanding spatial and temporal variabilities of severe hail and physical factors that contribute to those variabilities is important not only for advancing the science of extreme weather and climate, but also for mitigating hail damages. Past studies have shown that the spatial and temporal distribution of severe hail occurrence varies regionally and seasonally across the United States (Schaefer et al. 2004; Doswell et al. 2005; Cintineo et al. 2012; Allen and Tippett 2015). Changnon (1999) found that the severe hail frequency increased in the central Great Plains and southeastern United States and decreased in the Midwest and western United States over 1901–97. The increase in the central United States during the warm season was also seen during a 100-yr period, 1896–1995 (Changnon and Changnon 2000). These results revealed the importance of assessing the hail characteristics based on regions. In the U.S. southern Great Plains (SGP), most past studies have focused on the long-term hail trend and paid less attention to interannual variability.
Hail spatial and temporal variabilities are impacted by their associated large-scale environment, since large-scale circulation patterns influence the ambient conditions that favor severe storms. The subseasonal variability of hail occurrence in the contiguous United States was associated with the phase of the leading mode of the Madden–Julian oscillation (MJO; Barrett and Henley 2015) and the global wind oscillation (Gensini and Allen 2018). El Niño–Southern Oscillation (ENSO) has been found to modulate the hail and tornado occurrences in the United States during the winter (December–February) and spring (March–May) (Allen et al. 2015a, 2018; Cook et al. 2017; Childs et al. 2018). It has been shown that the probability of hail occurrence is larger in the cold phase of ENSO (La Niña) compared with the warm phase (El Niño), mainly through altering the large-scale environment. Besides ENSO, recent studies have also found that warmer sea surface temperatures (SSTs) in the Gulf of Mexico (GoM) correspond to more hail and tornado occurrences, and colder SSTs to fewer occurrences, during spring over the southern United States, mainly through influencing low-level moisture advection from the GoM (Molina et al. 2016, 2018; Molina and Allen 2019). The SSTs in the GoM are correlated with the ENSO neutral phase (Molina et al. 2016), but they are influenced by mesoscale features of the GoM such as the Loop Current and warm-core rings (Molinari 1987).
Aerosol loading may be another possible factor impacting severe storm activity in the central United States. Previous case studies have argued that biomass burning aerosols from Central America may contribute to the increase of hail and tornado occurrences over the central United States (Wang et al. 2009; Saide et al. 2015, 2016). Aerosols can impact convective intensity and hail formation and growth through aerosol–cloud interactions (i.e., acting as cloud condensation nuclei and ice nuclei; Rosenfeld et al. 2008; Khain et al. 2011; Loftus and Cotton 2014; Lebo 2014; Sheffield et al. 2015, Fan et al. 2018) and aerosol–radiation interactions (i.e., absorbing solar radiation and changing atmospheric temperature and humidity profiles and stability; (i.e., Tao et al. 2012; Fan et al. 2015). However, the significance of aerosol impacts on convective clouds has been found to depend on meteorological conditions such as convective available potential energy (CAPE), wind shear, and relative humidity (Fan et al. 2007, 2009; Khain et al. 2009; Storer et al. 2010; Lebo et al. 2012; Storer and Van den Heever 2013). In addition, aerosol concentrations can covary with meteorological conditions (Varble 2018) and variability in meteorological conditions may exert a stronger effect on convective properties relative to the aerosol effect (Grabowski 2018; Lebo 2018).
Climatological studies of hail have been mainly conducted using two types of data: 1) hail reports from National Oceanic and Atmospheric Administration (NOAA) Storm Prediction Center (SPC) and 2) radar-retrieved maximum expected size of hail (MESH). Compared with MESH, the primary benefit of hail reports is the greater confidence in hailstones occurring at the ground. However, there are spatiotemporal limitations and biases in hailstone measurement for the reports (Allen and Tippett 2015). Such limitations are influenced heavily by the population density (e.g., urban and rural regions), with more hail reports in densely populated areas compared to nonpopulated regions. Hail reports could also underestimate hailstone size (Blair et al. 2011), for example, due to surface melting prior to and during hailstone measurement (Blair et al. 2017). In the most recent decade or so, with the guidance of more weather radar observations, the hail report data can have an improved accuracy through verification of severe storm occurrence (Allen and Tippett 2015).
MESH estimates hail occurrence and size by using ground-based radar observations from the national Weather Surveillance Radar-1988 Doppler (WSR-88D) network (Crum and Alberty 1993; Crum et al. 1993). The radar network has nearly complete coverage over the continental United States (Maddox et al. 2002), and the data have been proven valuable for severe weather detection and warning decision-making. Compared to the hail reports, MESH provides more reliable and coherent spatiotemporal coverage (Lukach et al. 2017). However, due to the uncertainties relating radar observables to hail characteristics, it can have large biases in terms of the hailstone sizes (Ortega et al. 2009; Ortega 2018). Very recently, the algorithms for MESH were improved by Murillo and Homeyer (2019) and applied to the Gridded NEXRAD WSR-88D Radar (GridRad; Bowman and Homeyer 2017) dataset during the time period of 2004–16 in the central United States.
Given that the hail report data are the most robust for the most recent years with the guidance of dual-polarization radar observations, and the improved MESH data are also available, we analyze both datasets to understand the temporal and spatial variabilities of severe hail occurrences over the SGP region during the time period of 2004–16. These data are also combined with environmental and aerosol data to reveal the possible factors contributing to the observed severe hail variabilities. In addition, the cross-validation of both hail datasets is valuable for the severe weather community.
2. Data and methods
a. Hail datasets and statistical analysis methods
MESH was initially developed using an empirically driven power-law relationship fit between the 75th percentile hail size of 147 hail observations and the severe hail index—a measure based on the vertical integral of radar reflectivity (Witt et al. 1998). This study used newly developed MESH relationships to improve hail size discrimination. In particular, we used a power-law relationship fit to the 95th percentile hail size from approximately 6000 hail reports, which had better agreements with hail reports than a power-law fit to the 75th percentile hail size (Murillo and Homeyer 2019). The underlying data used to calculate MESH (GridRad) are on a regular longitude–latitude Cartesian grid with a horizontal resolution of 0.02° × 0.02° (about 2 km × 2 km), 1-km vertical resolution, and an hourly temporal resolution. Single-polarization and dual-polarization radar observations were combined from the Next Generation Weather Radar (NEXRAD) level II radar moments (volume scans) for the retrievals (Homeyer and Bowman 2017). The transition from single- to dual-polarization observations over the study period should not impact the MESH fields because the GridRad data are merged in a consistent way across all time periods and the dual-polarimetric variables are not used either 1) to construct the GridRad data or 2) to compute the MESH values. The transition from coarser resolution of the native WSR-88D polar grid (1 degree in azimuth by 1 km in range) to so-called “super resolution” (0.5° in azimuth by 0.25 km in range) in 2008 can impact MESH retrievals in regions with poor radar sampling, but such areas are largely excluded from analysis in this study. In particular, GridRad data include the number of radar scans that observed a grid volume at each analysis time and we only retain MESH data if at least 40 scans observed a given grid column.
Both hail datasets are used to analyze the temporal distribution of the total number of observed hail events over the SGP during spring (March–May) for the period 2004–16. For the annual trend of spring hail frequency, a 1° × 1° longitude–latitude gridded dataset was constructed using both the hail reports and MESH. The hail trend at each grid point was calculated using Theil–Sen slope estimator (Theil 1950; Sen 1968) and then smoothed by applying a two-dimensional 1.5-σ Gaussian kernel filter (Doswell et al. 2005; Allen and Tippett 2015) to reduce noises. Statistical significance of Theil–Sen slopes was assessed at 5% significance level using Kendall’s τ statistics.
For the statistical analysis of correlations of hail frequency with physical factors such as SST, ENSO, and aerosol loading, we applied both the least squares regression (LSR) with Pearson correlation coefficient and Theil–Sen regression (TSR) with Kendall rank correlation coefficient. The statistical significance was assessed at 5% significance level with two-tailed t test for LSR and with Kendall’s τ test for TSR throughout the paper. Since LSR analysis is sensitive to outliers but TSR is not, comparing LSR and TSR results can help us identify outliers.
b. Meteorological data
The analysis of synoptic environments is based on the North American Regional Reanalysis (NARR), which is available every 3 h at 32-km horizontal grid spacing and 45 vertical layers (Mesinger et al. 2006). The data at the closest time preceding each hail event are composited to determine the associated severe hailstorm environments. NARR adequately represents the convective environments, with realistic wind fields compared with observation (Gensini et al. 2014). We composited the following atmospheric variables: specific humidity, relative humidity (RH), temperature, and winds. The variables for convective environment such as CAPE, convective inhibition (CIN), and wind shear (0–6 km) are also examined. The maximum CAPE and CIN (i.e., based on surface air parcels) were obtained directly from the NARR dataset, and the 0–6-km wind shear was calculated by vertically interpolating winds at constant pressure levels to an above-ground level height coordinate (Gensini and Ashley 2011).
The global observed SST (°C) data used for exploring correlation with hail occurrences were obtained from NOAA’s Extended Reconstructed Sea Surface Temperature (ERSST) analysis version 5 (Huang et al. 2017). These data are produced on a fixed 2° × 2° longitude–latitude grid. To remove long-term warming signals, ERSST data were detrended before calculating SST anomalies using the 30-yr period of 1987–2016.
To analyze the relationship between hail occurrences and ENSO, the oceanic Niño index (ONI) from NOAA’s Climate Prediction Center was used. ONI identifies ENSO events in the tropical Pacific Niño-3.4 region (5°S–5°N, 170°–120°W). ONI is defined as the 3-month running mean SST anomaly relative to a 30-yr average roughly centered on each year (L’Heureux et al. 2013). The events are stratified into La Niña and El Niño phases based on the ONI threshold of ±0.5°C. We assess the statistical significance (95% confidence intervals) of the correlation between hail occurrences and ENSO using 10 000 bootstrapped resamples.
c. Satellite retrieved aerosol data and the approach to correlate with hail data
The quality assurance (QA) weighting average of Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua Level 3 aerosol optical depth (AOD) product (MYD08_D3, collection 6) at 1° × 1° longitude–latitude resolution is used in the analysis of correlation between hail occurrence and AOD, a proxy for aerosol loading. Terra and Aqua satellites cross the equator at around 1030 and 1330 local time each day, respectively. The MODIS data from Aqua are more appropriate for the analysis because Aqua overpass time is closer to hailstorm events since the majority of observed hail events over the SGP take place in the late afternoon and early evening (1500–2100 local time based on our analysis).
We processed the data as follows for the regression analysis between the hail data and AOD. First, we regridded the hail data from both hail reports and MESH to the resolution of the AOD data at each hour. Second, only AOD data at the same grids as the hail occurrence but prior to the hail occurrence time were chosen for the analysis. Since each AOD data point has the time of around 1330 local time, only the corresponding hail occurrences in the afternoon on the same day (before 0000 local time on the next day) are used. Last, the mean AOD in each year was weighted by the number of hail occurrences at each 1° grid, that is,
where n is the number of 1° grids containing hail events in a year and Nhail is the number of hail events in a 1° grid. This is to get an averaged AOD value over many 1° grids for each year with more weighting for the higher hail occurrence grids. With the weighting, the correlation between AOD and hail occurrence is slightly improved. We also excluded the lowest 5th percent of the significant severe hail days that are defined with the existence of significant severe hail during the PM hours (after 1200 and before 0000 local time on the next day) to avoid the small sample size problem (Doswell 2007). Note that the correlation over a large 1° grid would to some extent mitigate the effect of spatial shifts between aerosols and hailstorms.
We understand the uncertainties and caveats associated with using MODIS AOD data, but it is the only aerosol dataset available for such a long time period and over a large region, and it has been used for numerous observational studies in studying aerosol effects (e.g., Koren et al. 2010; Peng et al. 2016; Jiang et al. 2018; Zhao et al. 2019). Statistical significance of the correlation between AOD and hail occurrence is determined from 10 000 bootstrapped resamples.
d. HYSPLIT back trajectory model
To determine the origin of air parcels collocated with hail events and most likely ingested by hailstorms, we use the NOAA Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model, which is a complete system for calculating simple air parcel trajectories and is available on the NOAA Air Resources Laboratory website (http://www.arl.noaa.gov/). HYSPLIT is used to compute 48-h back trajectories initialized at altitudes of 1.0 and 1.5 km, since the surface air mass over the SGP region is mainly transported by low-level southerly winds below 1.5 km (Wang et al. 2018). Here we only show the results for the 1.5-km height on the hail days because the results for the 1.0-km height are similar. We calculate trajectories for the first hail event on every hail day since most of the backward trajectories for hail events on a given day followed a similar path. Winds used to drive the HYSPLIT model were obtained from the NCEP Global Data Assimilation System (GDAS) at a 1° longitude–latitude resolution.
a. Temporal and spatial hail variabilities and trends
We focus on the analysis for the significant severe hail category, which is defined as hail with a diameter greater than 2 in. (5.08 cm), since hail reports for such large hailstones are most robust (Blair et al. 2011; Allen and Tippett 2015). The results for the severe hail category, that is, 1 < diameter ≤ 2 in. (2.54 cm), are similar to those for significant severe hail.
Hailstones and tornadoes over the SGP are mainly associated with supercell thunderstorms (Blair et al. 2017) in the warm season (spring and summer). Our analysis confirms that hail (both significant severe hail and severe hail) occurs mostly from March to August and peaks in May in the SGP over 2004–16 based on both hail datasets (Fig. 1a). The majority of hail occurrences are in the spring season (March–May), which is consistent with previous results (Cintineo et al. 2012). In the summertime, the hail occurrences are generally low (Fig. 1a). Although June has moderate hail occurrences, SGP convective storms in summer occur within very different large-scale environments from those in the spring (Thompson et al. 2012; Tippett et al. 2014). In this study, our analysis focuses on spring hail events.
Interannual variability is substantial in spring, and there does not exist a statistically significant trend of spring hail frequency during the 2004–16 study period, as seen from both the hail reports and MESH datasets (Fig. 1b). We examine the temporal dependence of hail frequency using lag-one autocorrelation coefficients at the 5% significance level from two-tailed t tests. The lag-one autocorrelation coefficient is −0.04 for hail reports and 0.19 for MESH, suggesting a very weak temporal dependence of the hail frequency over the study period. Therefore, the temporally independent assumption used in our analysis is valid according to this statistical test. However, given the pronounced multidecadal variability of climate indices present in the United States and nearby regions, a 13-yr period is very likely not long enough to distinguish a greenhouse gas–induced long-term change from natural (internal) climate variability at the regional scale. The large interannual variability of hail frequency that we find is somewhat consistent with the study of Gensini and Ashley (2011) showing that climatologically the significant severe-weather environment varied both temporally and spatially in the SGP over 1980–2009.
Despite some inconsistencies in the spatial and temporal coverage of these two hail data sources, the interannual variability is significantly correlated between the two, with a Pearson correlation coefficient (r) of 0.81 and a null hypothesis (no correlation) p value of 0.001 based on a two-tailed t test. For MESH data, the number of invalid data (missing or bad radar data in the GridRad domain) has generally been decreasing during 2004–16 (Fig. 2a). To assess the impact of nonuniform sampling during the study period, we randomly subsampled the good data (spatially) at each time using the lowest sampling frequency that occurred during the study period (i.e., the minimum valid sample number among the 13 years) and repeated the interannual variability analysis. In that way, the number of samples analyzed from each year is effectively constant. The year with the minimum sampling number is 2004 (approximately 75% of the total grids, Fig. 2a). Surprisingly, this processing gives very similar results as the original results for both significant severe hail and severe hail (dark blue vs red in Figs. 2b,c). The resulting temporal distribution has a strong correlation with the dataset of all valid samples: r = 0.98 for significant severe hail (p = 0.000) and r = 0.99 for severe hail (p = 0.000). Thus, the differences of the radar sampling size among different years do not significantly affect the assessment of hail interannual variability and trend during the study time period. Therefore, we used the full dataset of MESH for the remaining analyses in this study.
By looking at the spatial distribution of the interannual trend for significant severe hail events, we see a statistically significant increasing trend in southern Oklahoma and northeastern Texas from both the hail reports and MESH based on the Theil–Sen slope analysis (Figs. 1c,d). The analysis of severe hail shows consistent results (Fig. 3). Changnon and Changnon (2000) found a similar increasing hail trend in Oklahoma and northeastern Texas, and a decreasing trend in southern Texas and southern Kansas during the 100-yr period of 1896–1995. However, this result is different from the proxy of severe tornado (EF2 or greater using the Fujita damage scale) annual trend revealed in Gensini and Brooks (2018), as an overall downward trend was found across the central and southern Great Plains over 1979–2017. The two hail datasets show inconsistent trends over Kansas, with an increasing trend from the hail reports and a decreasing trend from MESH, although these trends are not statistically significant. Allen and Tippett (2015) showed that the increasing hail reports over Kansas may be associated with the increasing reports from the observers such as storm chasers and research field projects for the period 1998–2014. Overall, both datasets indicate the spring hail frequency has large interannual variability over the SGP region, with a statistically significant increasing trend in southern Oklahoma and northeastern Texas during the past decade. About 9% of the grid points in the entire SGP domain have statistically significant trends at the 5% level using Kendall’s τ test for the hail reports (8% for MESH). We also applied a field significance test by using the resampling-based procedure (Wilks 2011), where we conduct the statistical test over the 13-yr study period using 1000 resampled replicates, to determine whether the fraction of grid points that show locally significant trends is greater than would be expected by chance (Livezey and Chen 1983). The field significance test indicates that the probability of finding 9% (and 8%) of the grid points in the SGP with statistical significance is only about 4% (and 3%) by chance, indicating that the field significance of the trends for the significant severe hail is statistically significant.
b. Factors contributing to hail interannual variability
1) Influence of SST anomalies over the GoM
Several previous studies have examined the influence of SSTs in the Pacific and Atlantic Oceans and GoM on hail occurrence in the central United States (Allen et al. 2015a; Barrett and Henley 2015; Molina et al. 2016; Lepore et al. 2017; Baggett et al. 2018). We follow those studies and examine SSTs in the GoM first using both least squares regression (LSR) and Theil–Sen regression (TSR). With LSR, the SST anomalies over the GoM appear to be moderately correlated with the interannual variability of hail events (r = 0.62 for the hail reports and r = 0.53 for MESH). We further divide the GoM into four subregions (G1–G4; Fig. 4a), and find that the SST anomalies in the northern GoM (G1 and G3) have the highest correlation (r = 0.73 for the hail reports and 0.61 for MESH). The goodness of fit, R2, for the univariate model from LSR is 0.54 for the hail reports (0.37 for MESH), meaning that the SST anomalies in the northern GoM might be responsible for 54% of the hail interannual variability observed from the hail reports (37% for MESH). With TSR, r = 0.54 for the hail reports and 0.46 for MESH, and R2 is 0.48 for the hail reports and 0.32 for MESH, showing a consistent result with LSR about the significant contribution of SST anomalies in the northern GoM to hail interannual variability (Figs. 4b–d).
The warmer (colder) SSTs over the northern GoM correspond to an increase (decrease) in hail occurrence over the SGP. By examining the relationship of hail occurrence over the SGP with the global SST variations over the same time period, we find a very weak correlation (r < 0.30, R2 < 0.1 using LSR), indicating that hail interannual variability over the SGP may be more influenced by the regional processes that alter SSTs over the northern GoM rather than the global environment. Those regional processes include, for example, the Loop Current (Sturges and Leben 2000; Ohlmann and Niiler 2005; Putrasahan et al. 2017). Molinari (1987) showed that the Loop Current, which varies seasonally and interannually, could influence SST variability by about 3°–4°C in the GoM.
To understand physical reasons leading to the observed relationship between the hail frequency variability and the SST anomaly in the northern GoM, we further analyze the synoptic environment differences between the four warmest (2006, 2011, 2012, and 2015) and four coldest (2007, 2010, 2013, and 2014) SST anomaly years using NARR. Figures 5a–d shows that the surface and low-level air is considerably warmer and moister over northern Mexico in the warmer SST years. This warmer and moister air is then advected into the SGP by the Great Plains low-level jet (LLJ; Fig. 5d). In addition, the LLJ is intensified by the enhanced pressure gradient from the Rocky Mountains toward the SGP (Figs. 6a,b), and the upper-level jet is stronger in the warmer SST years (Figs. 5d,e). The enhanced LLJ is coupled to the strengthened upper-level jet at 300 hPa through mass adjustments in the exit region of the jet stream (Uccellini and Johnson 1979) (Figs. 6c,d), resulting in larger wind shear (0–6 km) over the SGP (Fig. 5f), an important ingredient for severe hail-producing storms (Dennis and Kumjian 2017). Although the CAPE in most of the SGP region except southern Oklahoma and northeastern Texas is lower and convective inhibition is higher in warmer SST years (Figs. 5g,h), the significantly increased wind shear in combination with higher CIN suggests fewer, but potentially stronger convective storms (Trapp et al. 2007; Diffenbaugh et al. 2013).
2) Correlation with ENSO state
Besides SSTs, we also examine another large-scale environmental factor, which is the winter (December–February) ENSO state based on ONI, since the ENSO signal can be persistent from winter to spring in Niño-3.4 regions (Allen et al. 2015b; Lepore et al. 2017). We are not able to find a good correlation with seasonal mean ENSO phase with r = 0.10, R2 = 0.01 for the hail reports and r = 0.30, R2 = 0.09 for MESH using LSR when all of the years are considered (Figs. 7a,b). But with TSR, we find a better correlation with ONI (r = 0.23, R2 = −0.65 for the hail reports and r = 0.38, R2 = −0.31 for MESH; Figs. 7a,b), because TSR is insensitive to outliers, although the goodness of fit (i.e., R2) is worse since TSR is not aiming at minimizing sum squares of misfits. The different correlation values between LSR and TSR models and the comparably low R2 with both LSR and TSR models suggest the existence of one or more outliers. In fact, TSR results clearly identify 2016 to be an outlier (Figs. 7a,b). By removing the outlier 2016, the correlation between the spring hail frequency and ONI is much larger with r of 0.55 (0.56) and about 30% (31%) using LSR for the hail reports (MESH) (Figs. 7c,d). TSR shows r of 0.46 (and 0.52) and 20% (12%) with the statistical significance for the hail reports (MESH). Generally, the warm (El Niño) years as indicated by ONI > 0.5 correspond to low occurrences of the significant severe hail whereas the cold (La Niña) years (ONI < −0.5) correspond to high occurrences. Previous studies showed that during the La Niña phases, the upper-level jet stream and low-level moisture advection are enhanced compared with the El Niño phases, featuring a more favorable large-scale environment for convective storm development (Allen et al. 2015b; Cook et al. 2017).
3) Correlation with aerosol loading
To explore the potential contribution of aerosols to the observed hail interannual variability, we analyze AOD from the MODIS Aqua satellite and examine the relationship between the AOD before hail events and the interannual variability (Figs. 8a,b). When all 13 years are considered, both LSR and TSR show a certain level of correlation with AOD, with r = 0.52, R2 = 0.27 (r = 0.31, R2 = 0.14) for hail reports and r = 0.57, R2 = 0.32 (r = 0.33, R2 = 0.30) for MESH using LSR (TSR) (Figs. 8c,d). But the low R2 values suggest existences of outliers. The year 2016 was identified as an outlier in the ONI analysis and it is also identified as an outlier here. Indeed, AOD is low but hail frequency is high in 2016, which stands out as an obvious outlier. By excluding 2016, we find a stronger correlation of the hail frequency with AOD (Figs. 8e,f), with r = 0.76, p = 0.004 with LSR and r = 0.54, p = 0.016 with TSR for the hail reports (for MESH, r = 0.73, p = 0.007 with LSR and r = 0.55, p = 0.021 with TSR). The R2 value using LSR suggests that 57% of the interannual variability in the hail reports (54% in MESH) could be explained by the AOD differences. TSR shows very similar results, suggesting the robustness of the relationship between AOD and hail frequency. The R2 values are larger than those for the SST anomalies and ENSO, suggesting that the contribution of aerosols to the hail interannual variability might be greater. Worthy of mentioning is that we also examine the correlation between AOD from Terra and hail frequency, and the correlations (r = 0.48, R2 = 0.23 for hail reports using LSR without outlier) are indeed lower than that of Aqua (r = 0.52, R2 = 0.27), confirming that AOD from Aqua is more appropriate because the overpassing time is closer to hailstorm events in general.
Therefore, we find that GoM SSTs, ENSO, and AOD can be identified as factors associated with the interannual variability in spring significant severe hail occurrences. To examine their joint effects and respective contributions, we perform a multivariate linear regression with and without the outlier 2016 (Tables 1 and 2). The results show that the combination of the three factors can explain 63% (45%) of the hail interannual variability in the hail reports (MESH) using LSR with all 13-yr data considered (Table 1). By excluding the outlier 2016, the contribution rises to 75% (71%) for the hail reports (MESH), which is a lot higher than that accounting for the outlier (Table 2). Without the outlier 2016, both the hail reports and MESH consistently suggest that AOD may have the largest contribution (34% in the hail reports and 27% in MESH, respectively, using LSR), followed by SSTs over the northern GoM, and ENSO (Table 2).
While the apparent contribution of AOD variability to hail variability is high, aerosols can covary with meteorological conditions, which could contribute to the high correlation. To examine the covariability between AOD and meteorological parameters, we conduct a partial correlation analysis. The partial correlation measures the degree of association between two variables, with the effect of a set of controlling variables (meteorological parameters in this case) removed (Engström and Ekman 2010). If the partial correlations are similar to the corresponding total correlation, the meteorological parameters would have little effect on the total correlation. This technique was used to look at covariations of meteorological parameters with aerosols in a few previous studies (Zhao et al. 2018, 2019). The partial correlation calculations are conducted with nine meteorological parameters shown in Table 3 (top) eliminated individually and simultaneously. The differences between total and partial correlations are generally small (Table 3, top), indicating the observed correlation between the AOD and the hail interannual variability is not significantly related to the meteorological covariations. The largest difference is only 20% with all of the meteorological parameters removed for the MESH, suggesting that 80% of the correlation can be attributable to aerosol impacts. We further examine the correlation of AOD with the same sets of meteorological parameters as shown in Table 3 (bottom). The correlation coefficients are generally very low. Although r is 0.29 for U300 and 0.31 for CAPE in MESH, the p value is so large (0.52 for U300 and 0.33 for CAPE, indicating no statistical significance. Therefore, both the partial correlation analysis and the direct correlation of AOD with the meteorological parameters exclude a major role of the meteorological covariation in the observed correlation between AOD and the hail frequency.
Aerosols over the SGP in spring can be mainly transported from the regions to the south due to the predominantly southerly airflow (Wang et al. 2009). Previous studies documented the long-range transport of biomass burning aerosols from Central America to the south central and southeastern United States during spring (Wang et al. 2009; Saide et al. 2015, 2016). We examine the correlations of AOD in the SGP with that in Central America and find a relatively low correlation (r = 0.09, p = 0.781). The scatterplot shows a weak correlation and there is no significant outlier responsible for the low correlation (Fig. 9a). We further examine the correlations with AOD in Mexico, and do not find any good correlation (r = 0.29, p = 0.334; Figs. 9b, 10b). However, we find a high correlation with AOD in the northern Mexico (defined by the red contour in Fig. 10a) as shown in Figs. 9c and 10b (r = 0.62, p = 0.033). Since the aerosols we analyze for the hail events only, these results suggest aerosols from northern Mexico may contribute more to hail occurrences in spring over the SGP than those from Central America. Biomass burning and aged organic carbon are the two major aerosol types (comprised 40% and 31%, respectively) of the submicron aerosols in northern Mexico (Moffet et al. 2008). The low correlation with aerosols in Central America is likely related to the fact that a majority of low-level air (below 1.5 km) for the hail events cannot be traced to Central America based on the back-trajectory analysis. Our analysis is for hailstorms only, which does not exclude good correlations of aerosols in the SGP with those in Central America in other weather conditions.
4) Possible explanations for the outlier 2016
As mentioned previously, 2016 is an outlier in both AOD and ENSO contribution evaluations. It has a high hail frequency but a low AOD and a strong El Niño phase. The reason for low AOD in 2016 is because of the dominance of southwesterly and westerly low-level winds (Fig. 11a), bringing clean air (i.e., low AOD) into the SGP. The AOD in northern Mexico is actually high (Fig. 11b), but it does not influence the AOD over the SGP. With a strong El Niño phase in 2016, the synoptic environment should be less favorable for convective storm development in the SGP. However, we find that CAPE over the SGP in 2016 was higher and 0–6-km wind shear was stronger compared with those from the other three years (2008, 2011, and 2012) with similarly high hail occurrences (Figs. 12a,b).
To examine if the MJO is a potential factor leading to this more favorable meteorological environment in 2016, we conduct an analysis of MJO phase and corresponding meteorological conditions over the SGP. April and May 2016 are the focus of this analysis because they have the highest hail occurrences. The MJO phases are determined using the Real-time Multivariate MJO (RMM) index (Wheeler and Hendon 2004), and only the active MJO phases are examined. The active MJO phases are determined following Barrett and Gensini (2013). It is shown that, in April 2016, phase 8 is dominant and corresponds to the highest hail occurrences (Fig. 12c). However, based on Barrett and Henley (2015), under phase 8 of the MJO there is a very weak correlation (0.09) between hail-day anomalies and the synoptic environment, indicating that the large-scale environment is often unfavorable for hailstorm development. In May 2016, the dominant MJO phases are 1, 2, 3, and 4, and about 43% of significant severe hail events in spring occurring in those phases (Fig. 12d). In contrast, the corresponding dominant MJO phases are 5, 6, 7, and 8 in May of 2008 and 2011 and 8 in May 2012. We examine synoptic environment differences in the SGP region between the MJO phase 1–4 period in May 2016 and the mean of spring hail days of the other three high hail occurrence years. The CAPE over the SGP is smaller in 2016 and the 0–6-km wind shear is weaker, suggesting that MJO phases 1–4 do not correspond to the large-scale environment being more favorable for hailstorms than the average environment for high hail occurrence years as shown in Figs. 12e and 12f. Barrett and Gensini (2013) also found that tornado days were less likely in phases 2 and 3 compared with those in phases 5 and 8. In addition, Gensini et al. (2019) linked phases 8–1 of the MJO with the enhanced tornado activity in May 2019. However, the active phases of the MJO in 2016 are different from all these past studies. Therefore, our analysis suggests that the MJO may not have a strong influence on the high hail occurrences over the SGP in spring 2016.
The only positive correlation with high hail occurrences that we find is the relatively warmer SSTs in the northern GoM (Fig. 4b). However, the positive SST anomaly is still weaker than the other three high hail occurrence years. A possible contributing factor is abnormally high air temperature in 2016. The annual average temperature in 2016 for the contiguous United States was 1.2°C above the 1981–2010 average, which is the second warmest year since records began in 1985 (Blunden and Arndt 2017). However, more dedicated analyses are needed to develop a deeper understanding of factors contributing to the high hail occurrences observed in spring 2016 in the SGP.
For the increasing trend of spring severe and significant hail frequency observed in southern Oklahoma and northeastern Texas, we see increased CAPE and CIN, and stronger wind shear over these regions in warmer SST years, favoring stronger convective storms (Fig. 5). Since the SST warm anomaly years are dominant in the recent half decade (4 out of 6) while the SST cold anomaly years are dominant in the first half decade (Fig. 4b), the increase in the number of SST warm anomaly years is expected to be a factor contributing to the increasing hail trend at this region. Since the trend of AOD in southern Oklahoma and northeastern Texas is slightly decreasing (not shown), it is unlikely to be a factor contributing to the increasing trend in annual hail frequency.
4. Summary and discussion
With the most recent hail reports and improved radar-retrieved hail data, we have examined the occurrences of severe and significant severe hail over the SGP in a recent 13-yr period (2004–16) and the key factors that may contribute to the observed hail occurrences. This study focuses on the spring season (March–May; i.e., when hail occurs the most frequently). It is found that springtime hail occurrences in the SGP show a large interannual variability. Spatially, the increasing trend of hail occurrence in spring is statistically significant in southern Oklahoma and northeastern Texas.
We find three factors that appear to affect the interannual variability of hail in the SGP: aerosol loading, SST anomalies in the northern GoM, and ENSO. The interannual variabilities have a strong correspondence with SST anomalies over the northern Gulf of Mexico. 2016 is identified as an outlier for the correlations with both ENSO and aerosols loading. The three factors may explain ~63% of significant severe hail interannual variability in spring in the hail reports (45% in MESH) when all 13 years are considered, based on a multivariate linear LSR model. By removing the outlier 2016, the three factors can explain about 75% of the hail interannual variability for the hail reports and 57% for MESH. Without the outlier, both the univariate and multivariate linear models and both the hail reports and MESH datasets show that the correlation with AOD is larger than that with SSTs in the northern GoM and ENSO. Both the partial correlation and the direct correlation of AOD with the meteorological parameters show that the observed correlation between aerosols and hail interannual variability is not significantly controlled by the covariability with meteorological conditions. This indicates that aerosols may play an important role in the variability of significant severe hail occurrences over the SGP. Some possible mechanisms for increased hail production by aerosol loading can be increased supercooled water aloft that enhances riming by the increased cloud condensation nuclei (Khain et al. 2011), and enhanced convective intensity through larger latent heat release from condensation (Lebo 2014; Sheffield et al. 2015; Fan et al. 2018) as well as freezing processes (Rosenfeld et al. 2008). The effect of aerosols on hailstorms depends on aerosol properties and meteorological conditions such as cloud-base height (Carrió et al. 2014). Detailed process-level modeling studies are needed to study the significance of aerosol impacts and major mechanisms responsible for the effects.
We also find that the aerosol loading for the hail events generally correlates well with aerosol loading from northern Mexico because the majority of air parcels originate here as a result of prevalent southerly winds. Our analysis shows a weak correlation with aerosols in Central America because the majority of air masses near the hail events were not originating from or passing over Central America based on the back-trajectory analysis. Note our analysis here is for hail events only and does not exclude the possibility of a good correlation of aerosols in the SGP with those in Central America in other situations.
Warmer (colder) SSTs correspond to more (fewer) hail occurrences, because 1) warmer SSTs lead to higher moisture over the northern GoM and stronger 0–6-km wind shear over the SGP, and 2) the low-level jet is stronger by enhanced pressure gradient from the Rocky Mountains toward the SGP and stronger upper-level jet under warmer SST conditions, which enhances moisture transport into the SGP. This is consistent with several past studies that suggested warmer SSTs in the GoM are a significant source of moisture transport into the SGP, an important factor in influencing hail and tornadic storm development (Brooks et al. 2003; Jung and Kirtman 2016; Molina et al. 2016; Molina and Allen 2019). The increasing trend of annual hail frequency in southern Oklahoma and northeastern Texas may be the result of the increase in the number of SST warm anomaly years. ENSO can influence storm environments over the SGP through modulating the large-scale meteorological environment and SSTs over the GoM. Similarly, Lee et al. (2016) suggested a strong correlation of the ENSO phases with the large-scale processes conducive to tornadic storms in the United States.
Besides the hail occurrences, we also examined the relationships of hail days with ENSO, SSTs, and AOD to see how they are correlated. The correlations are all getting worse, with an r of 0.21, 0.15, and 0.12 for ENSO, SSTs, and AOD, respectively. We do not expect that there is a good correlation between hail days and aerosols, because hail days are largely determined by large-scale meteorology, which would not be affected by aerosols in the SGP. However, for a convective storm, the number of hail occurrences can be affected through aerosol impacts on cloud microphysics and dynamics.
We understand that observational studies can only reveal the relationship, and a high correlation does not mean causation. But such observational analysis is still a very useful first step, which can provide good insights for future work in modeling and field measurements. The potential relationships revealed here can be evaluated in greater detail by conducting a series of controlled model simulations, but that is beyond the scope of this study.
The study period in this work is only 13 years, which is likely too short to confidently distinguish an anthropogenic influence on the time series from the influence of possible natural multidecadal variability. Relatively short trends (e.g., 13 yr) over the United States may have large contributions from natural climate variability (e.g., McKinnon and Deser 2018; Wallace et al. 2015). It may also be too short to extract statistically robust correlations between hail and teleconnections (Brimelow et al. 2017). However, the results are obtained over the time period of the most reliable datasets that exist, which still may have significant implications regarding future hail occurrences. For example, Brimelow et al. (2017) suggested that future warming could lead to fewer small hail events but a more frequent occurrence of large hail over the SGP. The large variability in large hailstones and the high correlation with aerosols revealed in our study might add uncertainties as well as new perspectives in predicting the future changes of hail occurrences.
This study is supported by the U.S. Department of Energy Office of Science Early Career Award Program. PNNL is operated for the U.S. Department of Energy (DOE) by Battelle Memorial Institute under contract DE-AC05-76RL01830. C. Homeyer was supported by NSF under Grant AGS-1522910 and NASA under Award NNX15AV81G. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under contract DE-AC02-05CH11231.