A case study of South Australia’s severe thunderstorm and tornado outbreak 28 September 2016

2.1. Observational data

Satellite and radar observations, as well as ACCESS operational regional analyses (ACCESS-R; Puri et al. 2013; Bureau of Meteorology 2016b) are used to describe synoptic conditions associated with the severe thunderstorm and tornado outbreak of 28 September 2016 (Section 3) and for comparison with simulations (Section 4).

Radar images are derived from Adelaide’s Buckland Park S-band Dual Polarisation radar, at 0.5° elevation. Satellite data are processed by the Bureau of Meteorology from the geostationary satellite Himawari-8 operated by the Japan Meteorological Agency. Satellite images are also used to complement radar analysis for areas where storms and tornadoes were observed but were not well identified on radar due to poor or no radar coverage.

Locations of the seven reported tornadoes used in figures throughout this manuscript are based on the Bureau of Meteorology (2016a) report, with six of the tornadoes occurring during 06:00–06:30 hours UTC on 28 September 2016. In that report, tornadoes were confirmed based on damage assessment, damage reports and radar evidence.

2.2. Numerical model description

2.2.1. Set up of deterministic simulations

A case study of the 28 September 2016 event is conducted using a nesting configuration of the UK Met Office Unified Model (UM, ver. 10.6), which is the atmospheric component of ACCESS (Puri et al. 2013) used for research and operations in Australia. The model consists of a global model run (nominal 17-km horizontal grid spacing) that is nested down to grid spacing of 4.0 km, 1.5 km and 400 m (Fig. 1), with the inner domain size chosen to capture the area where seven tornadoes were reported. Each domain has 80 vertical levels, and the model top is 38.5 km.

Fig. 1. 

Outline of model domains (red), with horizontal grid spacing of 4 km for the larger domain, and 1.5 km and 400 m for the smaller domains.

All simulations use the Regional Atmosphere and Land Mid-latitude first release (hereafter RAL1-M) science configuration, developed for convection-allowing application of the UM (Bush et al. 2020). RAL1-M defines the science configuration of the physics and dynamics schemes of the atmosphere and uses Arakawa C-grid staggering in the horizontal (Arakawa and Lamb 1977) and Charney–Phillips staggering in the vertical (Charney and Phillips 1953). For a more detailed description see Bush et al. (2020).

The deterministic simulation is initialised at 15:00 hours UTC on 27 September 2016 (+9.5 h for local time) using initial conditions from global ACCESS (ACCESS-G), and all simulations are run for 24 h. Convection parametrisation is not used on the horizontal grid spacings of 1.5 km and lower as convection can be explicitly presented on these scales.

2.3. Model diagnostics

4.1. Overview of the 4- and 1.5-km simulations

Fig. 5 shows the upper-level (10-km) and near-surface (1-km) flow from the 4-km deterministic simulation from 6:00 hours UTC on 28 September 2016. Similar to Fig. 2, this shows a strong westerly jet extending over central and eastern South Australia and an upper-level low over the Great Australian Bight (Fig. 5a). At the surface a strong low-pressure system is situated over the Bight, with strong southerly winds to the south-south-west of the low centre. Associated with this low-pressure system is the north–south oriented cold front extending over central South Australia and collocated with a convergence zone with near-surface wind speed exceeding 20 m s⁻¹ (Fig. 5b). Analysis also shows that modelled MSLP at this time is 974 hPa, in good agreement with the observations (not shown).

Fig. 5. 

(a) Upper-level (10-km) and (b) near surface-level (1-km) winds (vectors and shading, m s⁻¹) and pressure (blue contours with a 2-hPa contour interval) at 06:00 hours UTC on 28 September (from the 4-km deterministic simulation).

To analyse how well the supercell thunderstorms and the QLCS are represented in the model, the 1.5-km deterministic simulation is analysed and compared with observations presented in Section 3. Fig. 6a shows simulated reflectivity around the time when the first thunderstorms associated with the outbreak began to develop along the cold front. Based on the Bureau of Meteorology (2016a) report and shown in Fig. 3a, the first thunderstorm cells formed on the northern Eyre Peninsula by 02:30 hours UTC on 28 September 2016 (12:00 hours local time). In the simulations, the first thunderstorm cells formed on the northern Eyre Peninsula at 02:30 hours UTC on 28 September 2016 (black box in Fig. 6a), in good agreement with the timing and location of the first observed thunderstorms. These thunderstorm cells continued to develop and formed a line while moving south-east ahead of the frontal boundary (Fig. 6b). Over the next 1.5 h (Fig. 7), the line of thunderstorms moved further east, extending in an approximately north-west–south-east direction over central South Australia.

Fig. 6. 

Composite reflectivity from the 1.5-km deterministic simulation at (a) 02:30 hours UTC and (b) 04:30 hours on UTC 28 September 2016. Black box in (a) denotes the location of the first thunderstorm cells that formed over the northern Eyre Peninsula.

Fig. 7. 

Composite reflectivity from the 1.5-km deterministic simulation at (a) 05:30 hours UTC and (b) 06:00 hours UTC on 28 September 2016.

This analysis shows some differences between the simulations and observations. For example, individual supercell thunderstorms associated with the reported tornadoes and seen in the satellite and radar image (cf. Fig. 4) are not fully resolved by the model. Also, the convective cells that extend further south are stronger than those observed. However, the overall timing, orientation and location of severe thunderstorms are captured reasonably well by the model. As mentioned in the introduction, numerical models with 1.5-km grid length (i.e. convection-allowing models) can more realistically represent convective storms and convective processes, but they are not able to fully resolve individual supercells or tornado-like signatures (e.g. Bryan et al. 2003; Hanley et al. 2016). In the following section, the 400-m simulation is analysed to determine if simulated storms are better resolved and if mesocyclones can be captured with this resolution.

4.2. Analysis of the 400-m simulation

Fig. 8 shows simulated radar reflectivity at 06:00 hours UTC on 28 September 2016, which is around the time when some of the tornadoes were reported (e.g. Bureau of Meteorology 2016a). The overall timing and location of the convective system is in good agreement with the observations (cf. Fig. 4). The northern cells (black box in Fig. 8) are stronger and better defined than the southern part of the convective system, indicated by differences in the simulated reflectivity (Fig. 8).

Fig. 8. 

Simulated radar reflectivity (dBZ) at 2 km at 06:00 hours UTC on 28 September 2016 from the 400-m deterministic simulation. The magenta dots are the observed locations of the tornadoes and the black box denotes the area shown in Fig. 9.

A close-up view of the simulated reflectivity and vertical velocity (Fig. 9) in the northern part of the system (black box in Fig. 8) reveals a hook-echo feature and a curved updraft (black box in Fig. 9) that coincides with the location of one of the reported tornadoes (northernmost magenta dot in Fig. 9). The hook echo indicates the presence of a mesocyclone (Davies-Jones 2015) at this location. The reported tornado is estimated to have formed at approximately 06:15 hours UTC on 28 September 2016 (Bureau of Meteorology 2016a; Sgarbossa et al. 2018), which is about 15 min later than in the simulation. Fig. 9 also shows that at the location of another reported tornado (the southern magenta dot in Fig. 9) larger values of reflectivity and vertical velocity are present, but there is no clear indication of a hook echo or cyclonic rotation. Similarly, analysis of the simulated reflectivity and vertical velocity in the southern part of the convective system, where other tornadoes were reported (cf. Fig. 8), does not show hook echoes or other rotational features (not shown).

Fig. 9. 

As in Fig. 8, but for the area indicated by the black box in Fig. 8. The black contours show upward vertical velocity (6–22 m s⁻¹, contour interval 4 m s⁻¹). The black box denotes the location of the hook echo.

Although a tornado cannot be resolved with the 400-m horizontal resolution, this analysis shows that the mesocyclone and features that indicate the possibility of a tornado (e.g. hook echo) are well captured. At this resolution the model can be used to investigate diagnostics such as UH (described in Section 2.3) to identify rotation and tornado potential in simulated storms. Although UH has been widely used as a diagnostic field in the United States National Weather Service, its use in Australia is relatively novel and to our knowledge, apart from the work of Sgarbossa et al. (2018), has not been widely used in operations.

Different thresholds for UH have been used in the scientific literature; for example, Kain et al. (2008) considered a threshold of 50 m² s⁻² as suitable for mesocyclone detection, whereas more recent studies used varying UH thresholds (75, 150 and ≥200 m² s⁻²), since higher UH thresholds are more suitable for reduced grid spacing (e.g. Clark et al. 2012, 2013; Gallo et al. 2016; Sobash et al. 2016b; Snook et al. 2019). As mentioned in Section 2.3, positive (negative) values of UH represent cyclonic rotation in the northern (southern) hemisphere. Here, a threshold of −120 m² s⁻² was used, chosen subjectively based on experimentation with different UH thresholds (not shown) to better depict mesocyclone path lengths.

Fig. 10 shows hourly minimum UH for the period between 05:00 and 07:00 hours UTC on 28 September 2016, which coincides with the time period of the reported tornadoes (Bureau of Meteorology 2016a). It shows elongated swathes of hourly minimum UH below −120 m² s⁻² and reaching values below −400 m² s⁻² in close proximity to the observed tornado tracks (magenta dots in Fig. 10). A long and coherent, vortex-like swathe of hourly minimum UH below −500 m² s⁻² coincides with the location of the northernmost observed tornado, indicating the mesocyclone that was identified in simulated reflectivity and vertical velocity (cf. Fig. 9). Interestingly, swathes of hourly minimum UH approximately coincide with locations of other observed tornadoes, where analysis of simulated reflectivity and vertical velocity did not indicate rotation; however, they are not as coherent as the northernmost swath of UH. Coherent tracks of high UH values may indicate potential tornado pathlengths but more research is needed to determine appropriate UH thresholds that best predict tornado pathlengths (Clark et al. 2012, 2013).

Fig. 10. 

Hourly minimum updraft helicity (m² s⁻²) for the periods (a) 05:00–06:00 hours UTC and (b) 06:00–07:00 hours UTC on 28 September 2016 from the 400-m simulation. The magenta dots denote the approximate location of observed tornadoes.

This analysis confirms that the use of UH as a diagnostic field would provide useful guidance for identifying the potential for tornado formation in this case. However, recent studies (e.g. Clark et al. 2012; Gallo et al. 2016; Sobash et al. 2016a) have shown that using solely UH as a tornado proxy can lead to large areas of false alarms in regions not always favourable for tornadogenesis.

Another parameter worth investigating as a tornado proxy alongside UH is the OW parameter (described in Section 2.3). Fig. 11 shows the OW parameter layer averaged between 1 and 4 km, at 4-min intervals for 05:28–06:28 hours UTC on 28 September 2016. It shows a coherent track of positive OW parameter that coincides with the location of the northernmost reported tornado, thus clearly identifying the mesocyclone in this simulation. In other locations further south, the OW parameter does not have coherent tracks, as they are more spread out and noisy. Other thresholds were tested when analysing the OW parameter, with higher thresholds showing only the northernmost track (not shown). In Fig. 11, a minimum threshold of 0.3 s⁻² is used. We also analysed the OW parameter for the 2–5-km layer (like UH) using the same threshold as in Fig. 11 and this also shows a northernmost coherent track of positive OW, but it is a lot noisier in other regions (not shown). This analysis shows that the OW parameter identifies the rotation associated with the mesocyclone in these simulated storms and gives a clear indication of a potential tornado path. This suggests that it can be used as an additional diagnostic alongside UH to assess tornado potential and identify regions where tornadoes may occur, thus potentially reducing false alarms.

Fig. 11. 

The Okubo–Weiss parameter (OW; positive values only), calculated as a 1–4-km layer average for the time interval of 05:28–06:28 hours UTC on 28 September 2016 from the 400-m simulation. Values of the OW parameter are spatially smoothed. The magenta stars denote the approximate location of the observed tornadoes.

4.3. Comparison with ensemble simulations

The analysis presented in Sections 4.1 and 4.2 is based on the deterministic, high-resolution simulation that enables the study of the dynamics and properties of tornadic storms, as well as highlighting how diagnostics such as UH and the OW parameter could be used as severe storm surrogates. However, it does not provide information on the uncertainty associated with the timing, location and intensity of the tornadic storms or tornado pathlengths (e.g. Hanley et al. 2016; Snook et al. 2019), which ensemble simulations can provide. Although tornado-resolving ensemble simulations are not practical due to the extreme computational costs involved (e.g. Zhang et al. 2015; Sobash et al. 2016b; Snook et al. 2019), convection-allowing (i.e. grid lengths of the order of O (1 km)) ensemble simulations are of considerable value for forecasting severe storms and associated hazards. In this section, 4-km and 400-m ensemble simulations are analysed using the model setup described in Section 2.2.

The ensemble simulations showed some differences to the deterministic simulation described above. There is an ~0.5–1-h difference in the timing of the event, likely due to the use of different global initial conditions in the ensemble simulations. Fig. 12 shows the simulated column maximum radar reflectivity for each ensemble member from the 4-km simulation. A substantial variation in the mode of convection is immediately apparent. Although members 1 and 2 display a QLCS structure quite similar to the deterministic forecast, member 5 displays weaker QLCS structure, member 4 is weaker throughout and members 0 and 3 have linear structures to the south of the reported tornadoes (magenta dots in Fig. 12) and discrete cells of convection to the north of the reported tornadoes.

Fig. 12. 

Column maximum simulated radar reflectivity at 05:00 hours UTC on 28 September 2016 from the 4-km ensemble simulation. The magenta dots denote the approximate location of the observed tornadoes.

Dial et al. (2010) found that the organisation of convection into linear or discrete modes is determined by the relative orientation of deep-layer shear vectors and mean wind vectors with respect to the initiating synoptic boundaries, such as cold fronts. The linear convective mode is favoured when shear vectors are oriented parallel to the synoptic boundary, whereas the discrete mode is favoured for shear vectors perpendicular to the boundary.

To understand the differences in convective mode between ensemble members (cf. Fig. 12), the orientation of surface–6-km shear vectors relative to the frontal boundary (in this case a cold front) is analysed. For brevity, only ensemble members 2 and 3 are presented here as they are good examples of linear and discrete convective modes, respectively. The storm motion vector is calculated based on Bunkers et al. (2014) and the surface–6-km bulk shear is calculated as the vector difference between winds at 6 km and winds at the surface (in this case 10 m).

Convection over the northern Eyre Peninsula develops a few hours earlier in member 2 than in member 3 (Fig. 13a, b) and along the frontal boundary (blue line in Fig. 13c), with shear vectors parallel to the frontal boundary (black ellipse in Fig. 13c) supporting a more linear mode of convection (Dial et al. 2010). Over time, convection in member 2 starts to resemble the QLCS mode that was observed (not shown). Convection in member 3 develops ahead of the front (Fig. 13b), with the linear mode favoured where shear vectors are more parallel to the boundary (closer to the low, black ellipse in Fig. 13d) and with the discrete mode to the north of the tornado reports where vectors are perpendicular (red ellipse in Fig. 13d). The approximate position of the frontal boundary was determined from 1.5-m temperature and MSLP analysis (not shown). This analysis agrees with the conceptual model of the relationship between the convective mode and the orientation of the shear vectors relative to the frontal boundary (Dial et al. 2010; Thompson et al. 2012).

Fig. 13. 

(a, b) Column maximum simulated reflectivity (dBZ) and (c, d) surface–6-km storm motion vector (black) and surface–6-km bulk shear vector (red) at (a, c) 01:30 hours UTC and (b, d) 04:40 hours UTC on 28 September 2016 from the 4-km ensemble simulation. The magenta dots denote the approximate location of the observed tornadoes and black and red ellipses in (c, d) denote where shear vectors are parallel or perpendicular to the frontal boundary (blue line).

A snapshot of the simulated radar reflectivity from the 400-m simulation (Fig. 14) shows more detail in the convective organisation between the six ensemble members. Convection in members 0 and 3 consists of both discrete cells and linear structure; however, the cells in member 3 are quite distinct and stronger than in member 0. Members 1, 2 and 5 display a QLCS structure similar to the deterministic forecast (cf. Section 4.2) but convection in member 5 is weaker than in members 1 and 2. Convection in member 4 is weaker throughout and patchy in the area of the reported tornadoes.

Fig. 14. 

Column maximum simulated reflectivity (dBZ) at 05:00 hours UTC on 28 September 2016 from the 400-m ensemble simulations. The magenta dots denote the approximate location of the observed tornadoes.

The results from the convective diagnostic tools from these simulations are similarly wide ranging, although we show only the UH for reasons of space. Hourly minimum fields of UH for the six ensemble members are shown in Fig. 15, with members 1 and 2 again being the most like the deterministic simulation. The four individual convective cells in member 3 each possess a strong band of intense UH, suggesting that these are individual supercells, in contrast to the QLCS structure of members 1 and 2. The weaker UH in the remaining members is indicative of less favourable, but not impossible, conditions for tornadogenesis.

Fig. 15. 

Hourly minimum updraft helicity (m² s⁻²) for the period 05:00–06:00 hours UTC on 28 September 2016, for the six ensemble members from the 400-m simulation. The magenta dots denote the approximate location of the observed tornadoes.