A Radar-Based Hail Climatology of Australia

The radar data used in this study are obtained from Level 1b of the Operational Radar Dataset in the Australian Unified Radar Archive (AURA, Soderholm et al., 2022). Our analysis encompasses all available data between 2000–2022, and mandates that each radar must contain over five years of operational observations in order to be included. Certain radars have also been excluded due to data quality concerns, such as geographically remote radars that did not historically transmit volumetric data, radars with persistent beam blockages, and various airport radars that are unsuitable for meteorological applications. These restrictions narrow our analysis to 52 radars, as illustrated in Fig. 1. Among the 52 radars examined, only eight are polarimetric¹¹1The Bureau of Meteorology’s Radar and Observation Network Uplift program dictates that all newly installed radars are polarimetric. At the time of writing, 16 of the 66 current operational radars are polarimetric, with this number expected to rise in years to come. and fourteen are S-bands (wavelengths $\sim$10 cm); both of which have important implications for our study, and merit separate discussions in Sections 22.2 and 3, respectively. The Level 1b dataset contains quality controlled radar data in polar coordinates, and due to the single-polarization limitations for most of the archive, the horizontally polarized equivalent radar reflectivity factor ($Z$, hereafter reflectivity) is used exclusively for hail detection. A filtering procedure is applied to remove ground clutter contamination (Gabella and Notarpietro, 2002; Heistermann et al., 2013), and radar reflectivity is calibrated using a ground clutter monitoring and spaceborne radar comparison technique (refer to Louf and Protat, 2023, for additional details). The latter is particularly important for our purposes, as calibration errors are known to severely impact hail size estimates (e.g., Warren et al., 2020). Recent three-way calibration checks for the Australian archive (between ground-based, ship-based and spaceborne radars, Protat et al., 2022) have found that post-correction calibration errors fall within $\pm$0.5 dB, theoretically leading to acceptable MESH errors of less than 5% (Warren et al., 2020, their Appendix A.).

Special consideration must be given for data quality issues when using radar data for climatological studies, particularly when the climatological variable in question may be strongly affected by weak nonmeteorological echoes, such as chaff, anomalous propagation and sidelobe contamination (e.g., Fabry et al., 2017). These effects are unlikely to impact our results, given the relative echo strength and the elevated nature of the observations used to derive MESH (>3 km altitude). A possible exception may arise when spurious high-altitude, high-reflectivity echoes are detected due to ground clutter from elevated topography, especially during atmospheric conditions conducive to anomalous propagation. We contend that these cases occur infrequently in the Australian network due to the absence of prominent topographic features (highest peaks $\sim$1500 m, Fig. 1), and direct the reader to the geometric justification given in Section 3.5 of Nisi et al. (2018) for a highly mountainous region. As in previous studies (e.g., Cintineo et al., 2012), a manual quality control procedure was applied to all identified hailstorms in the archive, which eliminated 1342 individual hail detections deemed spurious by a trained observer. We have determined that these evidently erroneous radar returns (such as complete volume saturation of >100 dBZ reflectivity) are likely attributable to hardware or signal processing failures and predominantly occur in the earlier portions of the archive for older radars.

The final data processing consideration involves the amalgamation of overlapping reflectivity observations from individual radars. Previous studies have highlighted the importance of integrating information from multiple radars in order to fill gaps in radar scanning strategies and compensate for range dependencies (Stumpf et al., 2004; Ortega et al., 2006). However, recent research by Brook et al. (2022) showed that gridding methods used to merge data into mosaics can degrade data quality, primarily through over-smoothing in the case of weighted average methods (e.g., Langston et al., 2007). Nearest-neighbor and linear interpolation methods can also propagate observational noise and introduce first-order discontinuities into the analysis (e.g., Lakshmanan et al., 2006). Single-polarized hail detection methods are highly dependent on the vertical distribution of reflectivity columns, the reconstruction of which is considerably degraded by both of the aforementioned gridding artifacts. Consequently, we have opted to regrid reflectivity observations using the variational interpolation method outlined by Brook et al. (2022) for each radar separately in our study, with a grid spacing of 1 km horizontally and 500 m vertically. This approach allows for quality control of individual radars without removing or recalculating the resulting mosaic when erroneous data from one radar is removed and ensures that isolated radars with no overlapping coverage are processed in the same manner as those that overlap. However, it also requires an additional post-processing step to merge the resulting hail detections for overlapping radars, which is detailed in Section 22.2. We deemed the regridding process essential in order to capitalize on its superior characterization of deep convective plumes, despite the considerable computational resources required for applying the variational approach to all hail events in the archive. These methodological improvements ensure that the gridded reflectivities are closest to those at the native resolution of the radar, where the original MESH constants were defined (Witt et al., 1998).

The Maximum Expected Size of Hail parameter (MESH, Witt et al., 1998) is used for hail detection in this study. Although originally developed as an object-based, single-valued parameter, we adopt the common approach of calculating MESH on regular grids (e.g., Stumpf et al., 2004). This process begins with the calculation of the kinetic energy flux based on the semiempirical relationship first introduced by Waldvogel et al. (1978),

We address the discrete temporal sampling rate of weather radars by implementing an advection correction algorithm from the PySTEPS package (Pulkkinen et al., 2019). This method interpolates MESH values along storm tracks by shifting values at time $t$ forward and values at time $t+1$ backward, before linearly weighting the two fields according to their temporal proximity (Anagnostou and Krajewski, 1999). Advection velocities are calculated across the domain by applying the Lucas–Kanade optical flow method to the MESH fields. The procedure is applied over contiguous time periods with valid hail detections (volumes with $\geq 30$ mm MESH), and temporal interpolation is performed at minute intervals. These time periods are calculated such that hail detections separated by less than or equal to 30 minutes are grouped together. This threshold was selected to accommodate earlier years of the archive that had 10 minute (rather than five minute) volume sampling times, ensuring hail streaks³³3We follow the terminology outlined in Changnon (1970) by referring to a hail streak as “an area of continuous hail with temporal coherence, which is considered an entity of hail generated within a thunderstorm”. A combination of spatiotemporally proximal hail streaks may then be subjectively grouped to form a hail swath for an event. from a single storm are not artificially split in cases with missing volumes due to hardware or communication failures.

Our analysis then departs from conventional MESH studies by adopting an object-oriented approach, where individual hail streaks are created by contouring interpolated MESH grids at the 30-mm level. Each hail streak is defined by the 30-mm contour polygon, and an example hail swath composed of hail streaks from three overlapping radars is shown in Fig. 2a. We conservatively (relative to prior studies, e.g., Soderholm et al., 2017a; Warren et al., 2020; Cintineo et al., 2012) limit swaths to a radial range of 20–140 km in our analysis to mitigate an underestimation of the true hail frequency due the radar’s scanning geometry (cone of silence and beam broadening). This conservatism is justified for the Australian radar archive, given the varying beam widths (between 1–2^∘, refer to Table C1) and the dual-PRF artifacts present beyond 140 km. A preliminary analysis of the range dependency of hail frequencies (not shown here) revealed that hail frequencies were relatively stable between our chosen range thresholds and no obvious underestimation was observed for radars with larger beam widths. We direct the reader to Bunkers and Smith (2013) for a more comprehensive discussion of these effects.

Hail streaks from individual radars are then merged by taking the geometric union of spatiotemporally overlapping polygons (refer to Figs. 2a and 2b), yielding spatially uniform, merged hail swaths. Merged hail streaks are assigned a range of statistics, including their area, start and end times, and optical flow velocities. Most of these statistics require an object-based approach, as they may not be easily calculated from an accumulated, gridded field (Nisi et al., 2018). The resulting merged hail archive (e.g., Fig. 2c), along with the associated object-based statistics, form the basis of our national and regional hail climatologies.

The precise areal definition of hail streaks in our object-based approach permits a novel, resolution-invariant method for calculating yearly or monthly hail frequencies. In a traditional approach using temporally accumulated (e.g., daily) MESH grids, hail frequencies can be evaluated at the original grid resolution by counting the number of times each pixel exceeds a mesh threshold (e.g., Cintineo et al., 2012). However, by chance, some isolated grid points may not experience a single hail storm throughout the archive duration, leading to the spurious claim that the true climatological frequency at that point is zero. One may attempt to account for the finite archive duration by applying various ad hoc smoothing techniques (e.g., Cintineo et al., 2012; Murillo et al., 2021), or calculating hail frequencies on a coarser spatial grid. Hail frequencies calculated using the latter approach (e.g., on a 10 km grid in Warren et al., 2020) yield higher values, as the number of hail detections within a 100 km² area is predictably greater than the number at each individual grid point. This resolution bias in hail frequencies prohibits reliable comparisons between studies with different grid sizes [e.g., Soderholm et al. (2017a) and Warren et al. (2020)], and confuses the true meaning of hail frequency estimates.

In this study, we define hail frequency as the number of times hail is expected at any single point within a grid box, per normalization period (e.g., yearly or monthly). We justify this choice as more intuitive to an observer, as it attempts to answer the common question: “How many times per year are you likely to observe hail?” This is distinct from: “How many times per year will it hail within 10 km from here?” This spatially invariant hail frequency is calculated by intersecting all hail streaks with a spatial region defined at any resolution. The cumulative area of all hail streaks within that region is then normalized by the area of the region ($A$) and its radar coverage duration ($N$, with quality controlled days removed). Formally, hail frequencies ($f_{h}$) are defined for each region as follows,

The irregular shape of overlapping radar coverage domains (see, for instance, Fig. 2a for three radars) further complicates the specification of hail frequencies in a regular grid. In grid boxes with overlaps in coverage, it is difficult to prescribe an $N$ value, as different regions within the grid box have been observed for different periods of time. This is especially true for larger grid boxes, such as the 100$\times$100 km grids illustrated in Fig. 2d. The solution is to intersect the radar coverage geometries with the underlying climatology grid, producing irregular “coverage regions”, each with a unique coverage duration. Hail swaths are then intersected with the coverage regions (black lines in Figs. 2d-f), and an individual hail frequency is calculated for each coverage region according to Equation 3. The result of these region-based frequency calculations is illustrated in Fig. 2e. Finally, we aggregate all region-based frequencies within each grid cell by preferentially weighting those regions with larger areas, longer archive durations, and more overlapping radars ($R$), as these factors all reduce uncertainty in the underlying frequency estimates. Weighted averages are computed as follows,

The Australian radar network covers roughly 30% of the country’s land surface, the distribution of which is illustrated in Fig. 1. As a result, some form of interpolation is required to create a continuous hail climatology throughout the country. This interpolation problem is particularly challenging due to large data voids, local variability in hail frequencies, and data noise introduced by the finite sampling period of the archive. Variational inverse methods are well suited to such a problem, as certain “prior knowledge” may be formally introduced to regularize an initially ill-posed system (e.g., Barth et al., 2014). Here, we seek to establish the following assumptions to help constrain the problem at hand:

1. 1.

   Radar-derived hail frequencies are a noisy but representative sample of the true climatological values.

2. 2.

   The true climatological hail frequency varies smoothly in space and time at our chosen resolution.

3. 3.

   No local extrema lie within data voids, and gaps are sufficiently constrained by the surrounding data.

4. 4.

   Hail frequencies reduce to zero at some distance from the coast.

We will briefly discuss the reasoning behind each of these assumptions before outlining their mathematical implementation. First, we justify the representativeness of our frequency estimates in Assumption 1 by noting the relative length of our archive (average $\sim$17 years, Appendix B) compared with seminal radar-based climatologies [e.g., four years in Cintineo et al. (2012) and eight years in Warren et al. (2020)]. We choose suitably sparse grid spacing to satisfy assumption 2 (100 km), which theoretically dampens high frequencies in the input data through spatial averaging and permits the use of smoothing constraints to “spread” valid information into data voids. Prior studies of the hail climatology in Australia (eg Allen and Allen, 2016; Bedka et al., 2018; Prein and Holland, 2018), provide justification for assumptions 3 and 4, where hail hotspots are not expected outside the current radar coverage or significantly out to sea. The extent to which the available data constrains the analysis within data voids will be investigated more thoroughly below.

We choose to perform the interpolation iteratively at two spatial scales: first at 200 km resolution to produce a provisional field, then at the final 100 km resolution using the former as a background field ($\varphi_{b}$). This type of nested analysis is common in variational problems (e.g., Laroche and Zawadzki, 1994; Dahl et al., 2019), and is useful in our case as valid data spreads more readily at a coarser resolution, providing a useful background for the final analysis. The cost function used in both cases is as follows,

Although $J_{s}$ may partially remove speckle noise from the frequency data, Brook et al. (2022) found that this constraint alone is unsuitable for denoising due to unwanted by-products produced when heavily weighting this constraint (oversmoothing and extraneous inflections in data voids, Shapiro et al., 2010; Brook et al., 2022). Complementary inclusion of an additional anisotropic total variation denoising constraint (Rudin et al., 1992) was found to adequately denoise the data, while damping extraneous inflections within the data voids and retaining high amplitude peaks within the analysis fields (refer to Appendix B in Brook et al., 2022). The denoising constraint is as follows,

Finally, the background penalty helps constrain regions outside the data coverage by penalizing deviations from an imposed background field ($\varphi_{b}$). At the provisional resolution (200 km), the background field is uniformly zero and only applies to grid points over the ocean through the application of a weighting matrix that increases exponentially from 0 over land to 1 out to sea ($W_{b}=\exp{-r^{2}/d^{2}}$) where $d$ is the distance from the coast and $r$ = 400 km controls the rate of weight increase). At the final resolution, the provisional retrieval is used as the background field, and $W_{b}$ is then set to 1 for grid points over land with no valid data. Formally, the background constraint is implemented as,

The aforementioned data voids may lead some readers to question the extent to which the surrounding data actually constrains the interpolated values within. In an attempt to satisfy these concerns and justify assumption 3, we apply the “clever poor man’s error” technique, which is an efficient method for assessing data coverage in variational interpolation problems (Troupin et al., 2012; Beckers et al., 2014). The procedure takes place as follows: set all available data to unit values, reduce the weighting constants by a context-specific factor, and perform the analysis as normal. Due to the difficulties involved in prescribing an exact covariance function for our interpolation procedure, we assume Gaussian covariance and consequently reduce the weights in the second step by a factor of $\sqrt{2}$ (Beckers et al., 2014). The resulting field approximates the “error reduction factor” (i.e., the extent to which the available data reduces analysis errors below background error state), where values close to one indicate the problem is well constrained by the data.

We treat the provision of weighting constants in Equations 6–8 as a parameter tuning problem (e.g., Brook et al., 2022, 2023). Suitable weights are derived experimentally and ultimately determine the relative strength of each constraint. The chosen values for our interpolation problem are as follows: $\lambda_{h}=0.5$, $\lambda_{d}=0.01$, and $\lambda_{b}=2$. The split Bregman optimization algorithm is used to minimize Equation 5 due to its applicability for mixed $\ell_{1}$–$\ell_{2}$ norms (limited to 5 inner, and 10 outer iterations; Goldstein and Osher, 2009; Ravasi and Vasconcelos, 2020).

We begin our analysis by evaluating how well the available radar data constrains hail frequency estimates in Australia. Figure Fig. 4a illustrates the radar network filled with unit data following the “poor man’s error” formulation (refer to Section 22.4). Significant data voids are apparent throughout the country, particularly in central WA, northern SA and western NSW. These areas are devoid of radars due to their sparse population and the fact they largely encompass Australia’s deserts. To test whether our interpolation technique is capable of propagating valid information into these voids, the error reduction factor was calculated as outlined in Section 22.4 and shown in Fig. 4b. In general, this quantity is equal to one over land and gradually decreases to zero away from the coastline, where the analysis is dominated by the background constraint and the radar data have minimal impact. The influence of the background constraint also extends slightly inland in regions with no radar coverage along the coastline (e.g., northeast and southeast WA). Such regions are less constrained by the radar data and should be interpreted with this in mind. It is possible to adjust the variational tuning parameters to ensure error reduction factors are close to one everywhere; however, we found that such combinations are prone to oversmoothing. We justify our choice of parameters by noting that the poorly constrained regions are small and that the resulting interpolated fields are a faithful representation of the underlying radar data (e.g., Fig. 5). Error reduction factors being close to unity almost everywhere over land validates our assumption that data voids are adequately constrained at our chosen analysis resolution.

Although Fig. 4b demonstrates that there is sufficient data coverage across the continent, the strength of that data coverage is highly heterogeneous. To quantify the strength of the data in regions with adequate coverage (error reduction factor > 0.9), we recalculate the error reduction factor for each grid point, this time setting the data at that grid point equal to zero. If the artificial data is placed in a data void, the analysis will be strongly influenced by it and the result will be close to zero. Conversely, if the point is inserted into a data-rich region, the addition of extra data will have less influence, and the result will be closer to one. The resulting field is therefore a quantitative assessment of the strength of radar coverage, with higher values indicating stronger coverage.

Fig. 4c reveals that the radar network strongly constrains the analysis around the coast, especially on the east coast. As noted earlier, the interior of the country contains large data voids, and the radar coverage strength is much lower there as a result. Interestingly, the coverage strength is not directly proportional to the distance to the nearest measurement. For example, small radar data gaps in southeast QLD and northeast NSW (Fig. 4a) do not cause a significant drop in this quantity, as they are surrounded by valid data on all sides. In contrast, the coverage strength in western NSW drops immediately when crossing into that data void, as it is unbounded on its northern and western edges. Regions with weaker coverage strength should be assigned lower confidence, as the data is sourced from more distant measurements, and the analysis is more sensitive to measurement errors. The 0.25 coverage strength contour from Fig. 4c is hatched in Figs. 5b and 6 to illustrate the fact that both data coverage and coverage strength should be taken into account when interpreting the subsequent national climatology.

Having discussed the intricacies of the chosen interpolation procedure, we present the resulting interpolated hail climatology in Fig. 5. The most notable feature is the broad hail frequency maximum located on the east coast of Australia, roughly situated south of Mackay in southeast Queensland and north of Sydney in northeast New South Wales. The highest hail frequencies are observed close to the coastline south of Brisbane, proximal to the Great Dividing Range (topography illustrated in Fig. 1). These findings are consistent with previous studies in the Australian region (e.g., Allen and Allen, 2016; Bedka et al., 2018). The position and size of this hail-prone region is similar to the hail environment-filtered overshooting top climatology in Bedka et al. (2018), but extends slightly further inland than their study suggests. This inland propagation is partially influenced by a high number of observed events in southwestern NSW (near the SA/NSW/VIC border, Fig. 5a). Manual quality control affirms the legitimacy of these events, but the frequency disparity compared to the surrounding region suggests that this feature may be a climatological anomaly. Although the interpolation procedure smooths such features towards the “climatological mean”, the anomalous information still propagates into the data void in western NSW. The limited data strength in this region (Fig. 4c), limits our confidence in how far west this hail maximum extends. The recent installation of two radars (Hillston and Brewarrina) largely fills the western NSW radar void and will help alleviate these uncertainties as more data become available in years to come.

Radar-based hail frequencies drop off appreciably south of Sydney, with annual frequencies on the order of 0.05 year^-1 across Victoria and southern South Australia. Hailstorms in this region are known to be associated with the passage of frontal systems, which can initiate thunderstorms in otherwise marginal environments. Zhou et al. (2021) demonstrated that these low instability and moderate to high shear hail environments are dominant in the high latitudes, where cooler temperature profiles also increase the likelihood of hailstones reaching the surface. Previous remote sensing studies have suggested that these relatively shallow hailstorms may result in an underestimation of hail frequencies in southern regions such as Melbourne and Canberra (Bedka et al., 2018; Punge et al., 2018). While there are numerous hail reports in these areas (e.g., Bedka et al., 2018), the reported sizes are smaller on average, reflecting less hail growth due to reduced precipitable water and weaker instability (Allen and Allen, 2016). As a result, a greater proportion of hailstorms and hail reports will be correctly overlooked by our 30-mm MESH threshold. It is also possible that MESH is poorly calibrated for detecting large hail due to these structural differences and underestimates the true frequency of hail in southern Australia (having been empirically fit to hail size reports in Oklahoma and Florida, Witt et al., 1998). Future work should investigate the relative accuracy of MESH estimates against hail reports for different high-population centers around the country.

On the west coast, radar coverage is limited to the coastline in the north but extends further inland in the Wheatbelt region between Perth and Kalgoorlie. Both regions exhibit elevated hail frequencies ($>$0.1 yr^-1), the first along the coastline between Exmouth and Broome and the second in the Goldfields region around Kalgoorlie. Both signatures appear fragmented in the observed radar frequencies (Fig. 5a), with disparate frequencies detected at adjacent grid points. Once again, the interpolation procedure tempers isolated high-frequency points to estimate the climatological mean. The result is a broad region of elevated hail frequency from Broome to Exmouth and around Kalgoorlie. Notably, the latter does not extend to the coastline in the south near Perth, Australia’s fourth most populous city. The land–sea contrast is especially apparent on the north coast of WA, where hail frequencies over the ocean are nearly zero, consistent with climatologically high values of convective inhibition (CIN) observed offshore in the region (Riemann-Campe et al., 2009). Other minor hail frequency maxima are observed by radars in central Australia near Alice Springs, and in northwest Queensland. The former region coincides with a collection of hail reports in the Australian Severe Storms Archive (Allen and Allen, 2016), although it is not identified in previous studies involving other remote sensing instruments (e.g., Dowdy and Kuleshov, 2014; Bang and Cecil, 2019). Bedka et al. (2018) did identify a local hail maximum in the goldfields region near the WA/NT/SA border, however, this is not mirrored in the radar archive, with the Giles radar west of Alice Springs showing little hail activity. The other local maximum near Burketown is also observed in the lightning climatology record of Dowdy and Kuleshov (2014), and is probably associated with low-level convergence in the Queensland heat low (Sturman and Tapper, 1996) and the onshore passage of cloud lines from the Gulf of Carpentaria (Drosdowsky and Holland, 1987).

Tropical radars north of Broome and Cairns (latitudes $\sim$18 and 17^∘S, respectively) exhibit minimal hail activity in our analysis, which is noteworthy for two reasons: 1) these findings may represent the first climatological application of MESH at such latitudes, and 2) previous hail climatology methodologies tend to overestimate hail occurrence in the tropics. These methods include satellite brightness temperatures (e.g., Bang and Cecil, 2019), overshooting tops (e.g., Bedka et al., 2018) and hail/severe storm proxies (e.g., Prein and Holland, 2018). All such methods are sensitive to tall, overshooting clouds in high-instability environments, which occur commonly in the absence of severe surface hail in the tropics. Several complementary factors are likely responsible for the absence of hail in the region. Firstly, vigorous warm-rain processes below the freezing level are thought to simultaneously limit the amount of supercooled liquid water and increase the number of competing hail embryos in the hail growth zone (Cotton et al., 2010). Second, the lack of wind shear in the tropics (e.g., Allen and Allen, 2016; Zhou et al., 2021) means that hail-producing storms are comparatively shorter-lived, as outflows from mature cells disrupt their own inflows (May et al., 2002). These structural differences likely reduce the potential updraft residence time, which is critical for hail growth (e.g., Kumjian et al., 2021).

Both structural and microphysical limitations discussed here dictate that hail observed aloft is mostly small (1–2 cm, May et al., 2002), and likely melts before reaching the surface due to typical temperature profiles in the tropics (Foote, 1984). To the contrary, Bang and Cecil (2021) demonstrated that spaceborne Ku-band radar reflectivity profiles are mostly consistent across tropical and subtropical hailstorms identified using different passive microwave detection methods. We propose that this discrepancy arises due to complex scattering and attenuation at Ku band radar frequencies, which reduces the reflectivity response for larger hail sizes relative to S-band and limits the information content available for hail size discrimination (Ryzhkov et al., 2013; Cao et al., 2013). A notable finding in Section 3 is the substantial impact of these effects at C-band frequencies, implying even more pronounced consequences at Ku-band. Additionally, the results presented here indicate that unlike Ku-band radars, ground-based S-band and corrected C-band radars effectively distinguish larger hydrometeors likely to reach the ground without melting, thereby minimizing the overestimation of hail frequencies in tropical regions. Admittedly, this analysis is limited to the available radar coverage, which excludes most of the Kimberley coast region in northwest Australia (see Fig. 1)⁷⁷7An important exception is the Broome radar that covers the Dampier Peninsula on the western fringe of the Kimberleys., where other remotely sensed hail (e.g., Bang and Cecil, 2019) and severe storm climatologies (e.g., Dowdy and Kuleshov, 2014) notably highlight as a hotspot. However, based on the findings in Section 44.1, we expect that the available radar data adequately constrains hail frequency estimates across tropical Australia. The ground-radar-based assessment of tropical hail frequency shown here appears realistic when compared to Severe Storm Archive hail reports (Allen and Allen, 2016) and aligns quantitatively with a 20-year climatology compiled from press reports by the Queensland Regional Meteorological Office (Frisby and Sansom, 1967), with the 1-in-10 year (0.1 yr^-1) and 1-in-20 year (0.05 yr^-1) isolines terminating around Mackay and Cairns, respectively.

The final result presented in this section concerns the seasonality of the Australian hail climatology. Fig. 6 presents interpolated monthly hail frequencies, calculated using normalized months to ensure an even sample size for each month. Segmenting the historical archive in this manner naturally limits the data available for each interpolation, potentially challenging our assumption that the observed data are a sufficiently representative sample of the actual climatological hail frequency. These limitations can lead to “patchy” hail frequencies for months such as November and February (Figs. 6e and 6h, respectively), which appear to be affected by climatological anomalies that arise due to the limited duration of the archive. We proceed in discussing the seasonal climatology with these limitations in mind.

While some hail activity is observed on the east coast as early as September, the hail season begins in earnest in October. Early season events are mostly concentrated in southern Queensland and progressively move south as the season advances. Note that the gap observed between Brisbane and Sydney in December and January likely reflects the lack of radar observations in the region (Fig. 5a). By February, the East coast hail frequency maximum has shifted south to around Sydney, showing a steady poleward trend in the intervening months. This aligns with global satellite observations made by Zhou et al. (2021), who suggested that the poleward shift coincides with increases in temperature and moisture at higher latitudes throughout the summer. Locally, Allen and Karoly (2014) and Warren et al. (2020) showed how these effects lead to a poleward shift in instability as the Austral summer progresses, particularly along the east coast. Conversely, the weakening of lapse rates in the absence of favorable upper-level support results in an overall reduction in hail-producing instability/shear environments in Queensland as the summer progresses (Allen and Karoly, 2014; Warren et al., 2020). As a result, the seasonal cycle is comparatively broader and later in NSW (November–March), a finding that is also reflected in hail reports from greater Sydney (Schuster et al., 2005).

Outside the east coast maximum, hail frequencies increase markedly throughout the country in November. This timing coincides with the arrival of continental heat troughs and their associated moisture transport to much of inland Australia (Sturman and Tapper, 1996; Hung and Yanai, 2004). This moisture, along with suitable upper-level support during late spring or early summer, provides favorable conditions for hail activity in those months, reminiscent of the seasonal climatology of convection in the Great Plains of the United States (Fabry et al., 2017). Under this hypothesis, the collective reduction in hail activity in January–February is a result of a warmer upper atmosphere as the summer progresses. Allen and Karoly (2014) also cited the importance of inland troughs in mid to late summer for storm development across inland Western Australia (WA). Hailstorms on the west coast of the country remain subdued until later in the season around February, consistent with the overshooting top-based climatology of Bedka et al. (2018). It is also worth noting that the only major hailstorm to impact the southwest capital city of Perth occurred similarly late in the season (March 2010, Buckley et al., 2010). Hailstorms persist in southern WA later than most of the country, likely due to their association with strong frontal systems that are capable of triggering hailstorms well into Autumn in the absence of strong instability (Allen and Karoly, 2014).