Horizontal Vortex Tubes near a Simulated Tornado: Three-Dimensional Structure and Kinematics

: Supercell thunderstorms can produce a wide spectrum of vortical structures, ranging from midlevel mesocyclones to small-scale suction vortices within tornadoes. A less documented class of vortices are horizontally-oriented vortex tubes near and / or wrapping about tornadoes, that are observed either visually or in high-resolution Doppler radar data. In this study, an idealized numerical simulation of a tornadic supercell at 100 m grid spacing is used to analyze the three-dimensional (3D) structure and kinematics of horizontal vortices (HVs) that interact with a simulated tornado. Visualizations based on direct volume rendering aided by visual observations of HVs in a real tornado reveal the existence of a complex distribution of 3D vortex tubes surrounding the tornadic ﬂow throughout the simulation. A distinct class of HVs originates in two key regions at the surface: around the base of the tornado and in the rear-ﬂank downdraft (RFD) outﬂow and are believed to have been generated via surface friction in regions of strong horizontal near-surface wind. HVs around the tornado are produced in the tornado outer circulation and rise abruptly in its periphery, assuming a variety of complex shapes, while HVs to the south-southeast of the tornado, within the RFD outﬂow, ascend gradually in the updraft.


Introduction
Supercell thunderstorms produce a wide spectrum of vortical structures with length scales ranging from the large, midlevel mesocyclone (5-10 km in diameter; [1,2]) down to vortices on the order of only a few meters, such as suction vortices [3][4][5]. One such manifestation of strong vorticity in supercells occurs in the form of three-dimensional (3D), elongated vortex tubes that are typically observed near tornadoes or in the periphery of their parent low-level mesocyclones [6][7][8]. They can occasionally be visually observed as condensation tubes or severely slanted funnels, given high-enough relative humidity and/or intense-enough cyclostrophic pressure drop inside them [9,10]. These vortex tubes are anecdotally referred to as "horizontal vortices" since, in contrast to tornadoes that are defined as vertical vortices, the axis of rotation of these vortex tubes is oriented primarily parallel to the ground.
Given the small scale and transient nature of horizontal vortices (hereafter referred to as HVs), they are hard to observe in Doppler radar data or to resolve in high-resolution numerical simulations; most evidence for their existence relies on videos or photographs [7][8][9]. Among the few existing observations, Wurman and Kosiba [5] provide evidence for large HVs south-southeast of two large tornadoes sampled by the Doppler-on-Wheels (DOW) mobile radar. In both cases, the HVs are located outside of the tornadic circulation but their sense of rotation (inferred from radial velocity plan-position indicators at two levels; see their Figure 14) is different; in one case the inferred horizontal vorticity vector points to the southwest while in the other it points to the north-northeast. This suggests different mechanisms may control HV formation. Houser et al. [9] provided an analysis of an HV interacting with the violent Piedmont-El Reno tornado on 24 May 2011 based on Rapid-scan X-Polarimetric Doppler (RaXPol) data and videographic observations. The HV, which is collocated with a weak-reflectivity band in reflectivity data, has a horizontal vorticity vector orientation similar to the Canton tornado case from Wurman and Kosiba [5], which pointed to the northeast and originated in a rear-flank downdraft (RFD) internal momentum surge to the south-southeast of the tornado, close to the surface. The HV wraps around the intensifying tornado, ascending in its circulation. The authors suggest possible mechanisms for the formation of the HV consistent with the observed horizontal vorticity vector orientation, which include: (i) baroclinic production along a warm RFD surge behind the primary RFD gust front, (ii) frictional torques in outflow air also behind the primary RFD gust front and (iii) reorientation of vertical vorticity associated with the tornado into a horizontal axis. Regardless of the mechanism, horizontal stretching of the vortex tube into the intensifying tornado is responsible for strengthening the HV to the point that cyclostrophic pressure drop in its core caused condensation to form.
Orf et al. [10], using a numerical simulation employing 30-m isotropic grid spacing to investigate the 24 May 2011 El Reno supercell and tornado, were the first to describe simulated HVs similar to visual observations. In the simulation, vortices ascend as funnel clouds on the periphery of the simulated tornado, in a similar manner to that described in Houser et al. [9].
As discussed by Houser et al. [9], surface friction is a plausible mechanism that can produce strong near-surface horizontal vorticity in the RFD outflow. This mechanism was shown earlier by Schenkman et al. [11] to have a significant impact on the vorticity budget of a developing tornado in a real case simulation of the 8 May 2003 Oklahoma City (OKC) tornadic storm [12]. In that study, surface friction acts on outflow and inflow parcels to produce strong near-surface horizontal vorticity that is abruptly tilted and stretched to produce pre-tornadic vertical vorticity centers that eventually coalesce into a tornado. This mechanism was explored in detail by Roberts et al. [13] and Roberts and Xue [14] in an idealized, single-sounding simulation of the 3 May 1999 Bridge Creek-Moore tornado. In summary, the authors show that surface friction acting on storm-induced flow produces near-surface crosswise horizontal vorticity that can be exchanged into streamwise vorticity as the flow bends cyclonically when converging toward a developing tornado (i.e., the river bend effect, e.g., Davies-Jones et al. [15]). The vorticity-rich parcels can then be tilted abruptly and stretched to produce strong near-surface vertical vorticity.
In order to better understand the behavior of HVs near strong tornadoes, in this study, data from a 100-m horizontal grid spacing, single-sounding simulation of a tornadic supercell are used to investigate the kinematics and dynamics of HVs near a simulated strong tornado. More specifically, our main objective is to provide a description of the 3D structure, types and evolution of HVs using direct volume rendering (DVR), a visualization technique used to generate high-quality displays of 3D numerical data offered by the Visualization and Analysis Platform for Ocean, Atmosphere, and Solar Researchers (VAPOR; [16,17]) software. The analysis is then validated by comparing the volume rendered vortical flows with visual observations of HVs in a real violent tornado. Finally, the kinematics of the near-ground flow in which HVs originate is briefly evaluated in order to infer potential formation mechanisms. More quantitative analyses of the dynamics within the simulation, including vorticity and circulation budget analyses following parcel trajectories and evolving material circuits, are deferred to future studies.
The reminder of the paper is organized as follows: Section 2 presents the numerical setup and methodology used to obtain the tornado simulation. Section 3 presents an overview of the simulated storm and associated tornadoes. Section 4 describes the structure and kinematics of HVs surrounding a simulated tornado using DVR of vorticity magnitude as compared to visual observations of real HVs, as well as traditional two-dimensional analyses of the vorticity field. Section 5 provides a discussion and summarizes the conclusions.

Numerical Experiments
The Advanced Regional Prediction System (ARPS; [18][19][20]) is used to perform the simulation presented in this study. The model domain is 84 km × 84 km in the horizontal, and 16 km in the vertical, with a Rayleigh damping layer applied above 12 km above ground level (AGL). The bottom boundary is flat. The horizontal grid spacing is 100 m and the minimum vertical grid spacing is 2 m at the surface, which stretches to 200 m above 10 km AGL, comprising 83 levels in the vertical.
The minimum vertical grid spacing of 2 m places the lowest scalar grid level at 1 m AGL. The reason for using such a small vertical grid spacing near the ground is twofold: (i) to provide increased vertical resolution very close to the surface and (ii) to place the first scalar model level closer to the surface than what is typically used in convective storm simulation studies (e.g., z ≥ 10 m AGL). Placing the lowest scalar level very close to the model surface is an attempt to increase the sample size of parcel trajectories that do not fall below the lowest scalar level, thus, avoiding the need to rely on extrapolated kinematic quantities [21]. Tornado and HV vorticity budgets along Lagrangian trajectories in this simulation will be examined in a future study, as indicated earlier. The drawback of using very small near-surface grid spacing is that the time step used in the explicit vertical advection scheme must also be very small; otherwise, the linear stability condition can be easily violated, especially when tornadoes are present. Therefore, a large (small) time step of 0.075 s (0.05 s) is used.
The model supercell storm is triggered by a thermal bubble within a horizontally homogeneous environment defined by a single sounding. The sounding is obtained from a full-physics 3 km grid spacing Weather Research and Forecasting (WRF) model simulation nested within a 15 km grid covering the contiguous United States for the 27 April 2011 devastating tornado outbreak in Mississippi-Alabama [26]. Conventional and radar observations were assimilated on the 15 and 3 km grids, respectively, using ensemble Kalman filter. The model sounding is extracted approximately 40 km southeast of the predicted storm corresponding to the Tuscaloosa-Birmingham tornadic supercell [7,27] at the 2100 Universal Time Coordinated (UTC) analysis time ( Figure 1). The analysis time is 1 h before the observed and predicted storms struck Tuscaloosa (at around 2200 UTC). The thermodynamic profile is obtained from a grid point over the city of Tuscaloosa (32.9 • N, 85.6 • W) and the wind profile is averaged over a 0.2 • latitude-longitude box centered on that grid point. A comparison of the surface air and dew-point temperatures at the Tuscaloosa airport (28.0 • C and 21.0 • C, respectively) with the WRF profile (27.8 • C and 21.9 • C, respectively) revealed a good agreement between the observations and WRF analysis. The location and time of the model-extracted sounding match the spatial and temporal ranges for proximity sounding suggested by Potvin [28] and minimize contamination from the predicted storm. A constant wind speed (u = 11 m s −1 , v = 17 m s −1 ) is subtracted from the sounding to keep the simulated supercell quasi-stationary during its tornadic phase. A similar procedure of extracting a sounding from a three-dimensional real-data simulation was taken by, e.g., Dawson et al. [24], as observed soundings at appropriate time and location were not available. Environmental hodograph of storm-relative winds averaged over a 0.2° × 0.2° latitude-longitude box centered over Tuscaloosa. The green arrow denotes the ground-relative storm motion vector that has been subtracted from the wind profile. Some relevant convective parameters are shown below the figures. The sounding and convective parameters were produced using the Sounding and Hodograph Analysis and Research Program in Python (SHARPpy) software [29].
Early, coarse resolution experiments (not shown) had spurious convection forming in the inflow of the developing supercell that eventually contaminated the storm. To avoid this problem, a small temperature increment was added to the sounding in the 720-2900 m layer, with a maximum 1.8 K added to the top of atmospheric boundary layer AGL. This increment rapidly decreases following a fourth-degree polynomial function such that at the bottom (top) of the modified layer, the perturbation is 0.2 K (0.5 K). The depth of the layer and the magnitude of the temperature increment were chosen by trial and error and were found to be the best option to eliminate spurious inflow convection while sustaining a vigorous tornadic supercell storm in the simulation. This procedure reflects the often-present top-of-boundary-layer inversion, and its direct effect was to reduce the mixed-layer Convective Available Potential Energy (MLCAPE) from 3442 J kg −1 to 3424 J kg −1 and increase the magnitude of Convective Inhibition (MLCINH) from 0 J kg −1 to 19 J kg −1 , respectively. Convection is initiated with an ellipsoidal thermal bubble perturbation centered at y = 66 km and x = 16 km with a horizontal (vertical) radius of 10 km (1.5 km). The maximum potential temperature perturbation at the center of the bubble is 6 K. This larger than typically used bubble amplitude is needed to develop a sustained supercell because of the larger inhibition added to the sounding as well as the large low-level wind shear in the environment.
The bottom boundary condition is semi slip with a constant drag coefficient Cd of 2.8 × 10 −2 applied to the surface horizontal momentum stress terms in the subgrid-scale turbulence parameterization scheme [13,14]. The value of Cd employed corresponds to an equivalent roughness length of 9.16 cm. The straightforward inclusion of surface drag in idealized convective storm simulations has the undesired effect of constantly slowing down the initial near-ground wind profile [30], destroying the base-state kinematic environment. In order to avoid this problem, the large-scale The red and green lines represent the temperature and dew-point temperature profiles, respectively, both in • C. The black dashed line indicates a parcel that ascends pseudo adiabatically from the surface. (b) Environmental hodograph of storm-relative winds averaged over a 0.2 • × 0.2 • latitude-longitude box centered over Tuscaloosa. The green arrow denotes the ground-relative storm motion vector that has been subtracted from the wind profile. Some relevant convective parameters are shown below the figures. The sounding and convective parameters were produced using the Sounding and Hodograph Analysis and Research Program in Python (SHARPpy) software [29].
Early, coarse resolution experiments (not shown) had spurious convection forming in the inflow of the developing supercell that eventually contaminated the storm. To avoid this problem, a small temperature increment was added to the sounding in the 720-2900 m layer, with a maximum 1.8 K added to the top of atmospheric boundary layer AGL. This increment rapidly decreases following a fourth-degree polynomial function such that at the bottom (top) of the modified layer, the perturbation is 0.2 K (0.5 K). The depth of the layer and the magnitude of the temperature increment were chosen by trial and error and were found to be the best option to eliminate spurious inflow convection while sustaining a vigorous tornadic supercell storm in the simulation. This procedure reflects the often-present top-of-boundary-layer inversion, and its direct effect was to reduce the mixed-layer Convective Available Potential Energy (MLCAPE) from 3442 J kg −1 to 3424 J kg −1 and increase the magnitude of Convective Inhibition (MLCINH) from 0 J kg −1 to 19 J kg −1 , respectively. Convection is initiated with an ellipsoidal thermal bubble perturbation centered at y = 66 km and x = 16 km with a horizontal (vertical) radius of 10 km (1.5 km). The maximum potential temperature perturbation at the center of the bubble is 6 K. This larger than typically used bubble amplitude is needed to develop a sustained supercell because of the larger inhibition added to the sounding as well as the large low-level wind shear in the environment.
The bottom boundary condition is semi slip with a constant drag coefficient C d of 2.8 × 10 −2 applied to the surface horizontal momentum stress terms in the subgrid-scale turbulence parameterization scheme [13,14]. The value of C d employed corresponds to an equivalent roughness length of 9.16 cm. The straightforward inclusion of surface drag in idealized convective storm simulations has the undesired effect of constantly slowing down the initial near-ground wind profile [30], destroying the base-state kinematic environment. In order to avoid this problem, the large-scale balance (LSB) technique discussed in Dawson et al. [31] is used. The fundamental idea in this technique is to assume that the initial horizontal wind profile is in a three-way balance among the horizontal-pressure gradient, Coriolis and frictional forces. Since the latter two forces are known from the initial wind profile, a pseudo-pressure gradient force is computed from the model time tendency of momentum equations at the first time step (i.e., before any modification caused by surface drag) for each model column. The pseudo-pressure gradient force (calculated for each grid point) is then added to the right-hand side of the horizontal momentum equations at all times at every grid point. The result is that the initial base-state wind profile remains virtually unchanged throughout the model integration.
Model integration is carried out using mode splitting, with the leapfrog (forward-backward) scheme for the slow (fast acoustic) modes. The model is integrated until 14,270.175 s using the 0.075 s large step. However, at this time, a tornado with strong near-surface updrafts is underway, and the vertical advection stability limit is violated since the first model layer is very shallow. From this time on, the simulation is restarted and run until 4 h 20 min using a large time step of 0.025 s while the small time step size is kept at 0.05 s. This is possible because vertical acoustic wave propagation is treated implicitly in the ARPS so that the small time step size is limited by horizontal grid spacing only. The large time step size is on the other hand limited by the very small vertical grid spacing in the case of large vertical velocity. The acoustic wave models are integrated using forward-backward integration schemes with time step size ∆τ, starting from the past time level t -∆t and ending the future time level t + ∆t, in n number of "small" time steps. When n = 1, ∆τ = 2∆t, hence the 'small' time step size is twice as large as the 'large' time step size [32,33].

Simulated Storm Overview
The simulated storm is fully developed into a supercell (e.g., with a hook-like appendage in the rainwater field and well-defined forward-and RFD gust fronts; Figure 2a) after 1 h of model integration as it gradually moves northwestward within the domain. This northwestward motion persists (but slows down) during the first 3 h of simulation. During this period, the storm develops only mesocyclone-scale vertical vorticity of O(10 −2 s −1 ) at low levels ( Figure 2a-c), with multiple low-level mesocyclone cycles taking place. However, tornado-intensity vortices do not develop until later in the simulation, a behavior that differs from the observed supercells in the 27 April 2011 tornado outbreak which were prolific tornado producers during a significant fraction of their life span. Several reasons may explain the differences between the simulated and observed storms as well the delayed tornado development in the simulation, which include: (i) inherent errors in the initial sounding extracted from the WRF analysis, (ii) the idealized, horizontally homogeneous design of the simulation that excludes heterogeneous environmental conditions encountered by the observed storm (e.g., complex terrain, cell mergers; [7,27]) and (iii) thermodynamic modifications applied to the initial sounding to prevent spurious convection. It should be pointed out, however, that early coarse grid spacing simulations (not shown) displayed a similar tendency for the storm to delay the production of strong surface vertical vortices in the first 2-3 h of integration, suggesting the delayed production of strong vortices may also be linked to the larger vertical shear present in the base-state wind profile [34].
After 3 h 16 min, significant structural changes occur in the supercell. Rapid strengthening of the low-level updraft at the tip of the hook (not shown) occurs coincidently with a sharpening (i.e., increase in confluence) of a left-flank convergence boundary (LFCB, as in Beck and Weiss [35]; compare Figure 2c,d). The increase in the degree of organization of the LFCB results in the formation of meridionally-oriented "vertical vorticity rivers" [35,36], which first appear at the base of downdrafts northwest of the low-level circulation (Figure 2d). These features, along with enhanced vertical vorticity in an RFD internal boundary, are thought to be the primary storm-scale sources of vertical vorticity to the developing low-level rotation in other simulations in the literature [35][36][37][38]. At 11,888.775 s, the first tornado forms at the tip of the hook as seen in the vertical vorticity field ( After 3 h 16 min, significant structural changes occur in the supercell. Rapid strengthening of the low-level updraft at the tip of the hook (not shown) occurs coincidently with a sharpening (i.e., increase in confluence) of a left-flank convergence boundary (LFCB, as in Beck and Weiss [35]; compare Figure 2c,d). The increase in the degree of organization of the LFCB results in the formation of meridionally-oriented "vertical vorticity rivers" [35,36], which first appear at the base of downdrafts northwest of the low-level circulation (Figure 2d). These features, along with enhanced vertical vorticity in an RFD internal boundary, are thought to be the primary storm-scale sources of vertical vorticity to the developing low-level rotation in other simulations in the literature [35][36][37][38]. At 11,888.775 s, the first tornado forms at the tip of the hook as seen in the vertical vorticity field (Figures 2d, 3a, and 4a). In this study, a tornado is defined as a persistent, concentrated vortex with surface vertical vorticity exceeding 0.1 s −1 and wind speed greater than 29 m s −1 , i.e., the minimum wind threshold of an Enhanced Fujita (EF) 0 tornado. presence of HVs [7]. The condensation funnel of this first tornado never touches the ground, a result of the insufficient cyclostrophic pressure drop within the tornado vortex. As the tornado matures, it briefly attains a maximum EF3 intensity at 12,032.550 s (Figures 2e, 3b, and 4b), with peak ground-level wind speeds of 62.2 m s -1 at 10 m AGL and surface core width (roughly estimated based on the highly-asymmetric radius of maximum wind at 1 m AGL) of about 200 m. After this short period of intensification, however, the tornado maintains only EF1 wind speeds during most of its life cycle. Throughout its life cycle, the tornado exhibits significant northeastward tilt with height [39] and is stronger near the ground during the tornadogenesis and maintenance phases (Figure 3a,b). In fact, nearly all the violent tornadoes on the 27 April 2011 tornado outbreak displayed pronounced northeastward tilt with height in addition to the almost ubiquitous presence of HVs [7]. The condensation funnel of this first tornado never touches the ground, a result of the insufficient cyclostrophic pressure drop within the tornado vortex.
In addition to the abovementioned LFCB vorticity river, less prominent rivers exist northeast of the tornado (Figure 2e), possibly associated with forward-flank convergence boundaries (FFDBs; [35]). The FFCBs merges with the LFCB just north-northwest of the tornado and may contribute to the vertical vorticity budget of the tornado (a topic that will be addressed in a future study). In addition to the abovementioned LFCB vorticity river, less prominent rivers exist northeast of the tornado (Figure 2e), possibly associated with forward-flank convergence boundaries (FFDBs; [35]). The FFCBs merges with the LFCB just north-northwest of the tornado and may contribute to the vertical vorticity budget of the tornado (a topic that will be addressed in a future study).
After approximately 17 min on the ground, the first tornado becomes completely occluded and wrapped in rain, resulting in broadening of the hook as the dissipating tornado moves rearward (southwestward) relative to the parent supercell [40][41][42] (Figures 2f and 4c). The near-ground vortex becomes completely detached from its upper portion and gradually decays (Figure 3c). The remnant vortex persists as a shrinking area of vertical vorticity embedded in precipitation that eventually mixes out within the turbulent downdraft outflow in the rear part of the storm. Despite that, the hook persists as a large low-level circulation, with an anticyclonic flare to its southeast (e.g., [43,44] After approximately 17 min on the ground, the first tornado becomes completely occluded and wrapped in rain, resulting in broadening of the hook as the dissipating tornado moves rearward (southwestward) relative to the parent supercell [40][41][42] (Figures 2f and 4c). The near-ground vortex becomes completely detached from its upper portion and gradually decays (Figure 3c). The remnant vortex persists as a shrinking area of vertical vorticity embedded in precipitation that eventually mixes out within the turbulent downdraft outflow in the rear part of the storm. Despite that, the hook persists as a large low-level circulation, with an anticyclonic flare to its southeast (e.g., [43,44]) ( Figure  2f). As the occluding downdrafts gradually fade, new vertical vorticity rivers form at the FFDB and LFCB as well as smaller feeders at the edge of a downdraft surge northwest of the low-level mesocyclone. A second tornado then forms at 13,950.225 s [45] (Figures 2g, 3d, 3g, and 4d). The second tornado, unlike the first one, becomes large (400 m wide at the surface) and strong enough such that the pressure drop within its core is able to produce a condensation funnel that extends all the way to the ground during its most intense phase ( Figure 3h). As this tornado matures (Figure 2h), it reaches ground-relative wind speeds at 10 m AGL of 86.7 m s −1 , corresponding to high-end EF4 intensity. The vortex is also tilted to the northeast but is wider at low levels than the first tornado at its maintenance phase (compare Figure 3b,e). Note also the well-defined "divided mesocyclone" structure [46] in Figure 4e, with the development of an occlusion downdraft east of the tornado.
Similar to the first tornado, the occlusion process wraps precipitation around the tornado, which eventually also becomes completely embedded in rain (Figures 2i and 4f). The circulation also broadens and becomes asymmetric, resulting in detachment of tornado and its associated LFCB and FFCB vorticity rivers, which now surround the occluded low-level mesocyclone. The decaying As the occluding downdrafts gradually fade, new vertical vorticity rivers form at the FFDB and LFCB as well as smaller feeders at the edge of a downdraft surge northwest of the low-level mesocyclone. A second tornado then forms at 13,950.225 s [45] (Figure 2g, Figure 3d,g, and Figure 4d). The second tornado, unlike the first one, becomes large (400 m wide at the surface) and strong enough such that the pressure drop within its core is able to produce a condensation funnel that extends all the way to the ground during its most intense phase ( Figure 3h). As this tornado matures (Figure 2h), it reaches ground-relative wind speeds at 10 m AGL of 86.7 m s −1 , corresponding to high-end EF4 intensity. The vortex is also tilted to the northeast but is wider at low levels than the first tornado at its maintenance phase (compare Figure 3b,e). Note also the well-defined "divided mesocyclone" structure [46] in Figure 4e, with the development of an occlusion downdraft east of the tornado.
Similar to the first tornado, the occlusion process wraps precipitation around the tornado, which eventually also becomes completely embedded in rain (Figures 2i and 4f). The circulation also broadens and becomes asymmetric, resulting in detachment of tornado and its associated LFCB and FFCB vorticity rivers, which now surround the occluded low-level mesocyclone. The decaying circulation then moves southwestward into the heavy precipitation of the hook and dissipates. Unlike the first tornado, however, as the second tornado becomes detached from its upper portion, multiple surface vertical vorticity centers (Figures 2i and 3f,i) associated with updraft pockets (Figure 4f) persist in the decaying broad surface circulation. After the second tornado dissipates, the simulated supercell maintains a well-defined hook in the rainwater field but no additional tornado forms until the end of model integration at 4 h 20 min (not shown). A dynamical analysis of the storm-scale processes leading to the formation of the tornadoes, their maintenance, and different decaying phases will be addressed in a future study.

Evolution and Kinematics of HVs near a Tornado
Throughout the development of the tornadoes in the simulation, the presence of HVs surrounding the tornadic circulations is ubiquitous, with striking similarities with the HVs presented in Orf et al. [10]. Using DVR, we now investigate the evolution, types and different structures of HVs exclusively focusing on the second simulated tornado since it is stronger and associated with a greater variety of HVs.

3D vorticity Structure and Visual Observations
During the intensification and maintenance phases of the tornado, there is enhanced activity of HVs near its outer circulation. It is also during this phase that these structures are more prominent and well defined. This is shown in a DVR of the 3D vorticity magnitude field in Figure 5, which highlights regions of vorticity magnitude > 0.15 s −1 , prior to and around the time the tornado attained EF4 strength. A supplemental video file SupVidPprtVort depicts the evolution of the vorticity magnitude field from 14,100.075 s through 14,298.000 s, when the tornado was rapidly intensifying. The general appearance of the vorticity magnitude field reveals a complex distribution of elongated vortices surrounding the tornado and its parent low-level mesocyclone. A wide spectrum of length scales is evident that includes vortices of nearly the same dimensions of the tornado down to scales near the grid spacing resolvability limit. This observation is supported by videographic evidence (Figure 6), which indicates that the width of some HVs can be as small as a few meters. (Note that the width of condensation tubes should be smaller than the actual HV circulation, e.g., [47]). The broad spectrum of HV scales (as well as other small-scale vortices in Figure 5) is a reflection of the turbulent nature of the storm's cold pool where they reside primarily, as also seen in the visualizations of Orf et al. [10].
Atmosphere 2019, 10, 716 9 of 19 circulation then moves southwestward into the heavy precipitation of the hook and dissipates. Unlike the first tornado, however, as the second tornado becomes detached from its upper portion, multiple surface vertical vorticity centers (Figures 2i, 3f, and 3i) associated with updraft pockets (Figure 4f) persist in the decaying broad surface circulation. After the second tornado dissipates, the simulated supercell maintains a well-defined hook in the rainwater field but no additional tornado forms until the end of model integration at 4 h 20 min (not shown). A dynamical analysis of the storm-scale processes leading to the formation of the tornadoes, their maintenance, and different decaying phases will be addressed in a future study.

Evolution and Kinematics of HVs near a Tornado
Throughout the development of the tornadoes in the simulation, the presence of HVs surrounding the tornadic circulations is ubiquitous, with striking similarities with the HVs presented in Orf et al. [10]. Using DVR, we now investigate the evolution, types and different structures of HVs exclusively focusing on the second simulated tornado since it is stronger and associated with a greater variety of HVs.

3D vorticity Structure and Visual Observations
During the intensification and maintenance phases of the tornado, there is enhanced activity of HVs near its outer circulation. It is also during this phase that these structures are more prominent and well defined. This is shown in a DVR of the 3D vorticity magnitude field in Figure 5, which highlights regions of vorticity magnitude > 0.15 s −1 , prior to and around the time the tornado attained EF4 strength. A supplemental video file SupVidPprtVort depicts the evolution of the vorticity magnitude field from 14,100.075 s through 14,298.000 s, when the tornado was rapidly intensifying. The general appearance of the vorticity magnitude field reveals a complex distribution of elongated vortices surrounding the tornado and its parent low-level mesocyclone. A wide spectrum of length scales is evident that includes vortices of nearly the same dimensions of the tornado down to scales near the grid spacing resolvability limit. This observation is supported by videographic evidence (Figure 6), which indicates that the width of some HVs can be as small as a few meters. (Note that the width of condensation tubes should be smaller than the actual HV circulation, e.g., [47]). The broad spectrum of HV scales (as well as other small-scale vortices in Figure 5) is a reflection of the turbulent nature of the storm's cold pool where they reside primarily, as also seen in the visualizations of Orf et al. [10].  The behavior of the most prominent HV in the simulation, i.e., the one highlighted by the black arrow in Figure 5, is further compared with visual observations in Figure 7. In the real Tuscaloosa tornado, a large HV is tilted upward in the right flank of the tornado as it revolves around it and then attaches to the cloud base in its forward flank (Figure 7a). A vertical cross section along the y-axis at 14,100.075 s in the simulation (Figure 7b) reveals the corresponding HV already located in the front sector of the tornadic circulation with horizontal vorticity predominantly oriented along the xdirection. In that sector, a lowering in the cloud base can be seen collocated with the HV. Although HVs are crudely resolved in the simulation, condensation driven by cyclostrophic pressure drop is still able to form when the HV nears the cloud base because of the near-saturation air at that level. The simulated perturbation pressure (Figure 7c) and vorticity magnitude (Figure 7d) fields bear a striking resemblance with the visual observation in Figure 7a, despite the more evident horizontal orientation of the HV in the simulation fields. The region where the HV attaches to the tornado in the perturbation pressure field (Figure 7c; the evolution of the perturbation pressure field from 14,100.075 s through 14,298.000 s is shown in supplemental video SupVidPprtVort) denotes a region where the upward tilting of the HV is reduced or even slightly tilted downward ahead of the tornado in the vorticity magnitude field (Figure 7d). Such reduced upward (or downward) tilt appears to be due to weaker updraft or downward motion in the secondary vertical circulation ahead of the tornadic circulation where the HV resides which, in this case, is enhanced by the pronounced northeastward tilt of the tornado. Thus, as this particular HV rotates around the tornado, it is affected by a minimum in upward motion (or downdraft) ahead of the tornado (i.e., along the tornado motion vector) and stronger updrafts at its northwestern and southeastern tips (transverse to the tornado motion vector; see Figure 4e). This pattern of horizontal gradient of vertical motion (∇hw) may result in the U-shape of the HV later in its life cycle (Figures 5b-e and 6c). In order to assess the similarity of the simulated vorticity structure with real-world observations of HVs, Figure 6 presents visual observations from video frames of the 27 April 2011 Tuscaloosa-Birmingham EF4 tornado, when several HVs were orbiting the tornado simultaneously. In both Figure 6a,b, the storm-relative viewpoint of the observers is approximately similar to that shown in Figure 5, looking from the southeast through the RFD region, but turning north and northwest (Figure 6c) as the tornado moved northeast. In both observed and simulated tornadoes, the HVs revolve counter-clockwise in the tornadic outer wind field and tend to align azimuthally around (but outside) the tornado core, thus, forming ring-like structures encircling the tornado. Although in a different context, this behavior bears remarkable resemblance to the process by which elongated vortices in isotropic homogeneous turbulence, the so-called "worms" (e.g., [48]), interact with a large, sustained columnar vortex in the DVRs of Takahashi et al. [49]. In their study, the interaction of the sustained vortex and the worms is shown to induce disturbances in the larger vortex core flow that can affect its dynamics. Given the vorticity arrangement shown in Figure 5, it seems plausible to hypothesize that vortex-vortex interactions analogous to those shown by Takahashi et al. [49] may occur between tornadoes and surrounding HVs.
Both observed and simulated HVs have a tendency to form in preferred storm-relative regions. The majority of the vortices appear or become well defined (as large tubes in Figure 5 or condensation funnels in Figure 6) to the rear (south and southeast) and right (east) sides of the tornadoes as they ascend in the outer circulation updraft. In fact, many HVs first appear close the ground in an arc extending from rear through right side of the tornado. This is seen in the convoluted vorticity distribution surrounding the simulated tornado (Figure 5a-c) and in the large, ascending vortex attached to the right flank of the observed Tuscaloosa tornado condensation funnel, shown in the leftmost arrows in Figure 6a,b.
The evolution of HVs shown in Figure 5 and in SupVidPprtVort highlights the diversity of structures and shapes the HVs can assume when interacting with the tornado. Large HVs close to the tornado initially tend to maintain their tornado-relative position (bottommost red arrow in Figure 5e and leftmost red arrows in Figure 6a,b), then eventually become detached from the surface as they are tilted by the updraft and finally spiral and evolve into complex shapes (e.g., uppermost red arrow in Figure 5e). Smaller vortices near the tornado edge, on the other hand, are rapidly captured in the tornado outer circulation and spiral upward, either maintaining their horizontal orientation or becoming severely distorted into a variety of shapes. Two such examples of complex shapes include the coil spring and U-shaped vortices in the forward sector of the tornado in Figure 5b,c and Figure 6a,b, denoted by the orange and black arrows, respectively. Clearly, the strong vertical motion gradients in the leading edge of the tornado (especially in forward-tilted tornadoes, as in both observed and simulated tornadoes presented here) and sinking motion ahead it must be involved in the creation of these highly-distorted vortex shapes (this aspect is further discussed later). We also note that vortices located farther away from the tornado core tend to be moved around with less change in shape or orientation (e.g., the U-shaped vortex in Figures 5d-f and 6c). This suggests that the rate of distortion of HVs is a function of the 3D shear strain rate surrounding the tornado core. Therefore, the length scale of the HVs and their distance to the parent tornado are both key to determining how their shapes and structures evolve near the tornadic wind field.
The behavior of the most prominent HV in the simulation, i.e., the one highlighted by the black arrow in Figure 5, is further compared with visual observations in Figure 7. In the real Tuscaloosa tornado, a large HV is tilted upward in the right flank of the tornado as it revolves around it and then attaches to the cloud base in its forward flank (Figure 7a). A vertical cross section along the y-axis at 14,100.075 s in the simulation (Figure 7b) reveals the corresponding HV already located in the front sector of the tornadic circulation with horizontal vorticity predominantly oriented along the x-direction. In that sector, a lowering in the cloud base can be seen collocated with the HV. Although HVs are crudely resolved in the simulation, condensation driven by cyclostrophic pressure drop is still able to form when the HV nears the cloud base because of the near-saturation air at that level. The simulated perturbation pressure (Figure 7c) and vorticity magnitude (Figure 7d) fields bear a striking resemblance with the visual observation in Figure 7a, despite the more evident horizontal orientation of the HV in the simulation fields. The region where the HV attaches to the tornado in the perturbation pressure field (Figure 7c; the evolution of the perturbation pressure field from 14,100.075 s through 14,298.000 s is shown in supplemental video SupVidPprtVort) denotes a region where the upward tilting of the HV is reduced or even slightly tilted downward ahead of the tornado in the vorticity magnitude field (Figure 7d). Such reduced upward (or downward) tilt appears to be due to weaker updraft or downward motion in the secondary vertical circulation ahead of the tornadic circulation where the HV resides which, in this case, is enhanced by the pronounced northeastward tilt of the tornado. Thus, as this particular HV rotates around the tornado, it is affected by a minimum in upward motion (or downdraft) ahead of the tornado (i.e., along the tornado motion vector) and stronger updrafts at its northwestern and southeastern tips (transverse to the tornado motion vector; see Figure 4e). This pattern of horizontal gradient of vertical motion (∇ h w) may result in the U-shape of the HV later in its life cycle (Figures 5b-e and 6c).
The general evolution of the HV shown in Figure 7 also suggests that, once the HV reaches the regions of strong ∇ h w ahead of the tilted tornado, it aligns perpendicularly to ∇ h w and attaches to this region. Ultimately, the HV becomes part of the horizontal vorticity field associated with the tornado-relative ∇ h w. This does not occur for weak, small HVs: as previously stated, they are quickly distorted by the strong velocity gradients as they approach the tornado outer edge. Therefore, the visualizations indicate that strong HVs can form far outside the tornado core via processes, such as frictional or baroclinic torques, and eventually become embedded in regions of sharp tornado-relative ∇ h w and attendant horizontal vorticity field. The general evolution of the HV shown in Figure 7 also suggests that, once the HV reaches the regions of strong ∇hw ahead of the tilted tornado, it aligns perpendicularly to ∇hw and attaches to this region. Ultimately, the HV becomes part of the horizontal vorticity field associated with the tornadorelative ∇hw. This does not occur for weak, small HVs: as previously stated, they are quickly distorted by the strong velocity gradients as they approach the tornado outer edge. Therefore, the visualizations indicate that strong HVs can form far outside the tornado core via processes, such as frictional or baroclinic torques, and eventually become embedded in regions of sharp tornadorelative ∇hw and attendant horizontal vorticity field.

Near-Ground Flow Kinematics and Potential HV Formation Mechanisms
Having analyzed the structure and types of HVs presented in the simulation, a brief assessment of the near-ground wind field surrounding the tornado and its associated horizontal vorticity is now presented. Even though vorticity budget analysis are not carried out in this study, the near-ground kinematic patterns around the tornado can provide valuable information about potential HV formation mechanisms and help guide future dynamical analyses.
As discussed earlier, among the variety of vortices observed in the simulation, a group of HVs is known to form near the ground in two key regions: in a circle around the base of the tornado and in an arc extending from rear through the right flanks of the tornado. The vortices near the base of the tornado continuously form and rapidly rise just outside the tornado core flow, while the ones in the rear and right flanks form longer near-surface tubes before accelerating toward the right side of the tornado and revolving around it, as previously described. One of the potential mechanisms for the rapid generation of HVs is via the Leibovich and Stewartson [50] instability. (An in-depth

Near-Ground Flow Kinematics and Potential HV Formation Mechanisms
Having analyzed the structure and types of HVs presented in the simulation, a brief assessment of the near-ground wind field surrounding the tornado and its associated horizontal vorticity is now presented. Even though vorticity budget analysis are not carried out in this study, the near-ground kinematic patterns around the tornado can provide valuable information about potential HV formation mechanisms and help guide future dynamical analyses.
As discussed earlier, among the variety of vortices observed in the simulation, a group of HVs is known to form near the ground in two key regions: in a circle around the base of the tornado and in an arc extending from rear through the right flanks of the tornado. The vortices near the base of the tornado continuously form and rapidly rise just outside the tornado core flow, while the ones in the rear and right flanks form longer near-surface tubes before accelerating toward the right side of the tornado and revolving around it, as previously described. One of the potential mechanisms for the rapid generation of HVs is via the Leibovich and Stewartson [50] instability. (An in-depth discussion of this mechanism can be found in Nolan [51].) Regardless of the genesis region, what is special about this class of vortices is that they typically erupt from the near-surface pool of enhanced horizontal vorticity. This shallow pool of horizontal vorticity that exists close to the ground (dark purple at the ground in Figure 5) is a direct result of surface drag and the resulting large near-surface vertical shear [11,13,14]. The tornado acts as a pump pulling the near-surface horizontal vorticity field upward: the vortices near the base of the tornado are abruptly displaced and/or tilted upward by the updraft while the ones in the right and rear flanks rise more gently along slantwise paths.
The idea that surface drag can be responsible for the generation of HVs can be linked to the predominant sense of rotation displayed by the vortices. In our 3D visualizations, virtually all large HVs, either embedded in the tornado outer circulation or RFD outflow, rotate around a horizontal axis with sinking motion ahead of them and trailing rising motion, as they orbit the tornado in the same way the HV presented by Houser et al. [9]. An analysis of the videos presented in Figure 5 also shows that the observed vortices in the Tuscaloosa tornado (as well as other tornado events, e.g., [7][8][9]) consistently display the same behavior. This suggests that HVs near the surface have horizontal vorticity vectors that point to the left of the horizontal wind vector at large angles, implying considerable crosswise horizontal vorticity [11,13,14,52]. This can be seen in Figure 8a, which shows the horizontal vorticity field at 1 m and 158 m AGL. Immediately outside the tornado core, in predominantly rotational flow, near-surface horizontal vorticity vectors, consistent with the near-surface HVs, point to the left of the prevailing horizontal wind in which they are embedded, thus having a large crosswise component. Above the surface (i.e., farther from the lower boundary; Figure 8b) the vorticity vectors become more streamwise, suggesting that HVs that arise from the near-ground pool of enhanced horizontal vorticity tend to exchange their crosswise vorticity into the streamwise direction, as they spiral around and upward near the tornado. This observation is consistent with generation of crosswise horizontal vorticity via surface friction in strong, accelerating horizontal flow [9,11,13,52]. discussion of this mechanism can be found in Nolan [51].) Regardless of the genesis region, what is special about this class of vortices is that they typically erupt from the near-surface pool of enhanced horizontal vorticity. This shallow pool of horizontal vorticity that exists close to the ground (dark purple at the ground in Figure 5) is a direct result of surface drag and the resulting large near-surface vertical shear [11,13,14]. The tornado acts as a pump pulling the near-surface horizontal vorticity field upward: the vortices near the base of the tornado are abruptly displaced and/or tilted upward by the updraft while the ones in the right and rear flanks rise more gently along slantwise paths. The idea that surface drag can be responsible for the generation of HVs can be linked to the predominant sense of rotation displayed by the vortices. In our 3D visualizations, virtually all large HVs, either embedded in the tornado outer circulation or RFD outflow, rotate around a horizontal axis with sinking motion ahead of them and trailing rising motion, as they orbit the tornado in the same way the HV presented by Houser et al. [9]. An analysis of the videos presented in Figure 5 also shows that the observed vortices in the Tuscaloosa tornado (as well as other tornado events, e.g., [7][8][9]) consistently display the same behavior. This suggests that HVs near the surface have horizontal vorticity vectors that point to the left of the horizontal wind vector at large angles, implying considerable crosswise horizontal vorticity [11,13,14,52]. This can be seen in Figure 8a, which shows the horizontal vorticity field at 1 m and 158 m AGL. Immediately outside the tornado core, in predominantly rotational flow, near-surface horizontal vorticity vectors, consistent with the nearsurface HVs, point to the left of the prevailing horizontal wind in which they are embedded, thus having a large crosswise component. Above the surface (i.e., farther from the lower boundary; Figure  8b) the vorticity vectors become more streamwise, suggesting that HVs that arise from the near- Another possibility for the horizontal vorticity of near-surface HVs is generation by baroclinic torques [15,24,[34][35][36][37]43,45]. However, if the horizontal vorticity in the HVs were mainly created baroclinically at the leading edge of a negatively buoyant RFD, the generated horizontal vorticity would point to the right of the downdraft and cold pool flows, giving an opposite sense of rotation than observed [5,43]. Yet baroclinity can still yield the observed vorticity if it is produced at the leading edge of an RFD warm surge, as suggested by Houser et al. [9]. Also vorticity thus generated should not be confined to be very close to the ground surface. This aspect is an important difference between Orf et al. [10] and the present study since the former employed a free-slip lower boundary condition while we include surface drag to capture the effects of surface frictional generation of horizontal vorticity. Recent studies [11,13,52] have found that surface-friction-generated vorticity can be an important or dominant source of vorticity the feeds a tornado vortex near the ground. This process is likely the main source of vorticity of the HVs observed in this study. More quantitative analyses of the simulation data will be needed to answer this question with any degree of certainty.

Discussion and Concluding Remarks
An idealized numerical simulation of a tornadic supercell is produced with the ARPS model at a tornado-resolving resolution and visualized with VAPOR; the data are analyzed, with the aid of 3D visualizations, to qualitatively determine the 3D structure, types and evolution of HVs near simulated tornadoes. The simulation includes surface drag with a semi-slip lower boundary condition in order to capture the effects of surface friction on the production of near-ground horizontal vorticity. Furthermore, the first scalar model level is placed at 1 m AGL to augment vertical resolution close to the surface to be resolve near-surface vertical shear, and to facilitate future trajectory-based analyses. The simulated parent supercell remains nontornadic during the first 3 h of life cycle, producing only mesocyclone-scale low-level vertical vorticity. After 3 h 16 min, the supercell transitions to its tornadic mode, producing two significant tornadoes. The first tornado lasts for 17 min and briefly attains EF3 strength, while maintaining EF1 winds for most of its life cycle. The second tornado is shorter lived but stronger, with peak winds reaching a high-end EF4 threshold. Both tornadoes exhibit significant tilt to the northeast with height.
HVs exist near the simulated tornadoes throughout their entire life cycle but are more prominent during the strengthening and maintenance phases of the tornadoes, consistent with the case study of Houser et al. [9]. A comparison of volume-rendered 3D vorticity magnitude fields and visual observations shows that the most significant HVs appear near the surface in an arc extending from the south toward the east sides of the tornadic circulation, embedded in both tornado outer flow and RFD outflow behind the gust front. Larger HVs tend to maintain their shape and position relative to the tornado longer than their smaller-scale counterparts, which quickly become distorted and are advected around and upward outside the tornado. Moreover, HVs farther away from the tornado tend to keep their shape and revolve around the parent circulation more slowly. This may be a result of reduced 3D shear strain rates outside and farther away from the tornado core flow. Once HVs are captured in the near-tornado, accelerating flow, very strong stretching occurs, which can occasionally form condensation inside the vortex tube via cyclostrophic pressure deficits, as shown in Figures 6  and 7 and also discussed by Houser et al. [9].
A closer inspection of a strong simulated HV reveals that the vortex is advected cyclonically toward forward flank of the tornado at low levels. Because the tornado tilts forward with height, strong horizontal gradients of vertical motion exist ahead of it. Our visualizations suggest that strong HVs can form within the storm cold pool and later become attached to the zone of strong horizontal gradient of vertical motion ahead of a forward-tilted tornado. Thus, HVs can occasionally be assimilated in the horizontal vorticity field associated with tornado-relative vertical motion gradients. In fact, many of the fast moving, forward-tilted tornadoes observed on the 27 April 2011 tornado outbreak displayed HVs in their forward sectors during a significant portion of their life cycle [7].
Another relevant aspect seen in both DVRs and videos of the Tuscaloosa tornado is the tendency for HVs to have horizontal vorticity vectors with a large crosswise component at low levels, similar to the 24 May 2011 El Reno tornado studied by Houser et al. [9]. Regions of large horizontal near-ground vorticity collocated with and pointing to the left of strong near-surface horizontal winds strongly suggest a possible role by surface drag in the production of HVs. Recent studies have shown that surface friction can produce significant near-ground horizontal vorticity within RFD outflows [11] and in the storm's inflow [13,14] that is ultimately tilted and stretched to produce tornadic vortices. It may well be possible that the HVs described in this study are coherent-structure manifestations of the vorticity-generation processes discussed by Schenkman et al. [11] and Roberts et al. [13], revealed here owing to the use of 3D visualizations and comparisons with real tornado footage. Recently, 3D visualizations of simulated tornadic supercells have provided crucial insights into storm-and tornado-scale process involved in tornado formation and dynamics [10,53]. Mechanisms other than surface friction, such as baroclinity along RFD warm surges behind their primary gust fronts and reorientation of near-tornado vertical vorticity into the horizontal [9], may also be important for the generation of the variety of 3D vortex tubes in the vicinity of a tornadic wind field. The complex distribution of HVs near the simulated tornadoes in the RFD outflow suggest all these mechanisms may be acting together. A conceptual model of the 3D structure, relevant flow features and what we believe to be the most plausible mechanisms for HV formation is presented in Figure 9.
that surface friction can produce significant near-ground horizontal vorticity within RFD outflows [11] and in the storm's inflow [13,14] that is ultimately tilted and stretched to produce tornadic vortices. It may well be possible that the HVs described in this study are coherent-structure manifestations of the vorticity-generation processes discussed by Schenkman et al. [11] and Roberts et al. [13], revealed here owing to the use of 3D visualizations and comparisons with real tornado footage. Recently, 3D visualizations of simulated tornadic supercells have provided crucial insights into storm-and tornado-scale process involved in tornado formation and dynamics [10,53]. Mechanisms other than surface friction, such as baroclinity along RFD warm surges behind their primary gust fronts and reorientation of near-tornado vertical vorticity into the horizontal [9], may also be important for the generation of the variety of 3D vortex tubes in the vicinity of a tornadic wind field. The complex distribution of HVs near the simulated tornadoes in the RFD outflow suggest all these mechanisms may be acting together. A conceptual model of the 3D structure, relevant flow features and what we believe to be the most plausible mechanisms for HV formation is presented in Figure 9.  The vertically-oriented vortex represents an intensifying or mature tornado while slantwise, detached vortex tubes represent more horizontally-oriented vortices in the periphery of the tornado. Regions of enhanced frictional generation of horizontal vorticity in strong, near-ground horizontal wind embedded in the RFD outflow are shaded in purple. Magenta arrows in the middle panel show the tilted circulation on the forward side of the tornado. Representative vortex lines associated with the HVs are displayed in light orange, with circular arrows indicating their sense of rotation. Strong surface RFD flow is indicated by the curved black arrows. Right column: The cloud field associated with a tornado producing HVs consistent with the vorticity field displayed in the left column and visual observations. The RFD-related clear slot is annotated. The storm is moving to the northeast, as indicated by the red arrow the upper panel.