Remarkable progress has been achieved in understanding the vorticity source responsible for tornadogenesis. Nevertheless, the answer to this question remains elusive, particularly after introducing surface friction in realistic tornado simulations. In this study, a simulation using the Weather Research and Forecasting (WRF) Model is conducted based on the F3 supercell tornado that hit Tsukuba City, Japan, on 6 May 2012. The simulation uses triply nested domains, and the tornado is successfully reproduced in the innermost domain with 50-m horizontal grid spacing. The circulation analyses reveal that the frictional term is the dominant vorticity source responsible for the vortices at both the pretornadic and tornadogenesis times. The detailed vorticity source analyses of the air parcels show that the vorticity of the tornado at the genesis time mainly originates from the frictionally generated crosswise vorticity near the ground. The crosswise vorticity is directly tilted (or first exchanged into streamwise vorticity and then tilted) into vertical vorticity when the air parcels enter the tornado. A rear-flank downdraft (RFD) surge from the south and west sides of a primary low-level mesocyclone (LMC) may trigger tornadogenesis by increasing the convergence near the ground. The RFD surge is not necessarily associated with the baroclinically generated vorticity. In this study, the baroclinity is weak across the hook echo, which may cause a lack of baroclinically generated vorticity in the RFD surge.
To understand the mechanism of supercell tornadogenesis, identifying the responsible vorticity sources is a vitally important issue. Despite several decades of study on this issue, a complete understanding remains ambiguous. Previous studies included field observations and idealized numerical simulations. It was suggested that the baroclinic effects are prominent in generating the vorticity of tornadoes (e.g., Davies-Jones and Brooks 1993; Adlerman et al. 1999; Straka et al. 2007). The baroclinic effects are due to the horizontal buoyancy gradient associated with the downdrafts in a storm. Horizontal vorticity is generated by baroclinity when airflow approaches a tornado along the rear-flank downdraft (RFD; Straka et al. 2007; Markowski et al. 2008), along the forward-flank gust front (FFGF; Wicker and Wilhelmson 1995; Markowski et al. 2012; Markowski and Richardson 2014), or at the periphery of downdrafts aloft (e.g., Dahl et al. 2014).
Although idealized simulations have traditionally neglected surface drag, several new studies in the past several years have included parameterized drag. Roberts et al. (2016) suggested that surface drag generates crosswise vorticity, and the crosswise vorticity is consequently converted into vertical vorticity in a tornado. Markowski (2016) showed that the contributions of environmental and storm-generated (baroclinic or frictional) vorticity to the circulation of a near-surface mesocyclone depend on different environmental wind profiles and supercell stages.
Realistic tornado simulations are becoming more frequent in recent years. Mashiko et al. (2009) reproduced a minisupercell tornado in a rainband of Typhoon Shanshan using the Japan Meteorological Agency Nonhydrostatic Model (JMANHM). Environmental streamwise vorticity was suggested to be the main component of the vorticity sources responsible for tornadogenesis. The environmental vorticity is associated with the strong low-level veering shear in the typhoon. Schenkman et al. (2014) utilized the Advanced Regional Prediction System (ARPS) to examine a supercell tornado that occurred in the 2003 Oklahoma City storm. The vorticity sources of four pretornadic vortices that led to two tornadoes were investigated. The crosswise vorticity generated by surface drag was found to be responsible for all the vortices. The crosswise vorticity was exchanged into streamwise vorticity and consequently tilted into vertical vorticity. Mashiko (2016b, hereafter M16) simulated the 2012 Tsukuba tornado and proposed that frictionally generated circulation is the dominant source of the pretornadic vortex, while baroclinically generated circulation is prominent for the tornado at genesis time.
In addition to the controversy regarding the vorticity sources of tornadic vortices in different studies, recent studies indicate that even for a given tornado, the vorticity sources can be different with different model settings. Yokota et al. (2018, hereafter Y18) simulated the same 2012 Tsukuba tornado by ensemble forecasts and analyzed the circulation sources of the vortices at the times of tornadogenesis and maximum intensity. The authors revealed that circulation generation is mainly (not always) due to the frictional terms. The baroclinic terms can also be dominant at times.
The RFD surge (also named the internal outflow surge or secondary RFD) behind the primary rear-flank gust front (RFGF) was found to be related to tornadogenesis in some recent studies. These studies suggested that the RFD surge enhances convergence, which leads to tornadogenesis (e.g., Kosiba et al. 2013; Mashiko et al. 2009). However, the generation of vorticity associated with the RFD surge is still under investigation. Marquis et al. (2012) suggested that the secondary RFD may baroclinically generate vorticity and contribute to the tornado at the maintenance stage by vortex tilting. M16 suggested that the baroclinically generated vorticity on the tip of the hook echo is transported into the tornado via the RFD surge. Schenkman et al. (2014) suggested that the frictionally generated vorticity is prominent in the RFD surge.
Previous studies have also found that the RFD surge can induce cyclic mesocyclogenesis (Adlerman et al. 1999; Adlerman 2003; Beck et al. 2006). The developing RFD wraps around the primary mesocyclone and causes the gust front to surge outward. This surging gust front initiates a new updraft on the downshear flank of the primary mesocyclone. Then, a new mesocyclone is formed in the new updraft.
The Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008) is a widely used meteorological model. However, the application of this model in resolving tornadoes has been rarely demonstrated. Most studies have been limited by coarse grid resolutions. Xu et al. (2015) simulated a bow echo system with 0.8-km horizontal grid spacing. The significant contribution of surface friction to the generation of low-level mesovortices was emphasized. Marquis et al. (2014, 2016) investigated the 5 June 2009 Goshen County tornadic supercell. Vorticity budget analyses were conducted on several sets of air parcels from a near-surface mesocyclone.1 The authors suggested that the baroclinity within the forward-flank region is the main source of vorticity of the near-surface mesocyclone throughout its lifetime. A more recent study (Flournoy and Coniglio 2019) simulated a tornadic mesovortex in a quasi-linear convective system with 333-m horizontal grid spacing. Three separate airstreams entered the mesovortex and contributed to the vertical vorticity in different ways. In the southerly inflow, the tilting of crosswise vorticity is predominant, while in the northerly inflow and a rear inflow jet, the tilting of streamwise vorticity is predominant. Sun et al. (2019) simulated an EF4 supercell tornado in Funing, China, and discussed the sensitivity of flow features in the tornado to grid resolutions. The authors recommended that horizontal grid spacing of 50 m or less is needed for studies of tornado dynamics.
As mentioned above, the discussion of the vorticity sources of a tornado remains controversial particularly after introducing frictional effects. The main purpose of this study is to understand the mechanism of tornadogenesis by considering the sources responsible for the vorticity. To achieve a detailed understanding of the tornadogenesis process, a realistic simulation by WRF with large-eddy simulation implementation (WRFLES; Moeng et al. 2007) is conducted on the 2012 Tsukuba tornado. The finest horizontal grid spacing is 50 m, and the tornado is explicitly resolved. The vorticity generation mechanism of the flow entering the tornado is explored by backward trajectory with a particular focus on an RFD surge. Y18 indicated the uncertainty in the tornadogenesis mechanism for a given tornado. Thus, the generated tornadoes in M16 and Y18 and the present study can all be regarded as solutions of the Tsukuba tornado. The atmospheric condition resolved in the present study is physically possible and favorable for tornadogenesis. In this context, the analysis of tornadogenesis mechanism in the present study is expected to achieve meaningful results.
Section 2 introduces the experimental design. Section 3 describes the general characteristics of the mesocyclone and the explicitly resolved tornado. Section 4 introduces the methods of circulation and vorticity analyses. Section 5 discusses the circulation sources of near-ground vortices at the pretornadic and tornadogenesis stages. Vorticity budgets along the trajectories of air parcels are investigated in section 6. A discussion and conclusions are presented in section 7.
2. Experimental design
The present study focuses on the Tsukuba tornado, which is an F3 tornado that occurred on the afternoon of 6 May 2012 (Japan Meteorological Agency 2012). The WRF Model (version 3.7.1) simulation has triply one-way nested domains (Fig. 1a). The outermost domain (WRFLES1km) covers most of the area of eastern Japan using 300 × 300 horizontal grid points with 1-km spacing. This domain aims to resolve the environment of the storm and runs from 0000 to 0600 UTC. Nested within WRFLES1km, the second domain (WRFLES250m) starts from 0200 UTC with 250-m horizontal grid spacing. The number of horizontal grid points of WRFLES250m is 801 × 701, and the domain covers most of the Kanto Plain, where the convective storm developed. Starting from 0230 UTC, the innermost domain (WRFLES50m) with 50-m horizontal grid spacing is nested within WRFLES250m. The number of horizontal grid points of WRFLES50m is 1501 × 1501, and the tornado is explicitly resolved (Fig. 1b). WRFLES1km and WRFLES250m have 62 vertical levels with spacings that stretch from approximately 40 m near the ground to approximately 640 m at the model top with a constant pressure of 60 hPa. WRFLES50m has 81 vertical levels in which the spacing near the ground is 10 m, and there are 24 levels below the 1-km height. A fifth-order horizontal advection scheme and a third-order vertical advection scheme are used. A third-order time-split Runge–Kutta scheme is used to advance the computation in time.
The model uses the Morrison et al. two-moment bulk microphysics scheme (Morrison et al. 2009). Cumulus parameterization is turned off. Radiation is handled through the RRTMG scheme (Iacono et al. 2008). The five-layer thermal diffusion land surface scheme and the revised MM5 Monin–Obukhov surface layer scheme (Jiménez et al. 2012) are also included. The Monin–Obukhov similarity theory (MOST) is based on steady and well-mixed flows. However, the outflow of a thunderstorm is normally unsteady, which would hinder the reliability of the surface layer scheme. Markowski et al. (2019) suggested that the MOST tends to underpredict the surface layer wind shear in the storm environment, especially in the outflow of the storm.
In recent years, the “gray zone” problem (Wyngaard 2004) has been encountered in downscaling simulations using mesoscale models. When the grid size is at an intermediate scale between the LES upper limit (slightly less than the integral turbulence scale) and the mesoscale lower limit (no resolved turbulence), neither the planetary boundary layer (PBL) schemes nor the LES models are appropriate. For the simulation of a supercell storm, Fiori et al. (2010) suggested that even for 1-km horizontal grid spacing, the LES closure produces better results than the PBL schemes and “gives acceptable representation of storm-scale supercell structure and evolution”. We also reach similar conclusions in sensitivity tests (not shown). Moreover, considering that the simulation period is daytime, and the atmospheric condition is unstable in the storm environment, it is reasonable to expect deep mixing. Thus, the present study uses an LES model. The LES model is based on the 1.5-order turbulent kinetic energy (TKE) closure scheme (Deardorff 1972). A drawback of the LES model should be noted here. Markowski and Bryan (2016) showed that when simulating a homogeneous, quasi-steady boundary layer using LES, unrealistically large wind shear would develop if the turbulence is not sufficiently resolved (especially in a laminar flow). The environmental inflow of a storm is likely error prone due to this phenomenon. In the present simulation, turbulence has partly developed in the environmental inflow, which may alleviate this problem. Another problem with the LES model is that the 1.5-order TKE model may overestimate the wind shear on the lowest model levels because the resolution of these grids is inadequate to resolve the turbulent eddies near the ground (Mason and Thomson 1992). Note that the overestimation of wind shear caused by the two problems may be partly offset by the underestimation caused by the MOST scheme.
The topography and land cover of WRFLES1km are interpolated from the U.S. Geological Survey (USGS) dataset, the resolution of which is 30 arc-s. For WRFLES250m and WRFLES50m, a 100-m meshed dataset provided by the Geospatial Information Authority of Japan is used. The JMA mesoscale analysis data are utilized for the initial and boundary conditions of WRFLES1km. These data are updated every 3 h and have a 5-km horizontal resolution. Data assimilation is not used in the present study.
3. Overview of the mesocyclone and tornado in the WRFLES simulation
Figure 2 shows the supercell generated in WRFLES250m. A midlevel mesocyclone (MMC; defined at the 3-km height) is associated with the hook echo of the supercell. The MMC has a diameter of approximately 3 km and a central pressure deficit of approximately 8 hPa. The storm translational velocity is (17, 6) m s−1 (zonal and meridional components), which is defined by the average moving velocity of the MMC center from 0258 to 0308 UTC. The MMC center is defined at the location of maximum vertical vorticity. Although the 250-m horizontal grid spacing simulations can hardly explicitly resolve small structures, such as tornadoes, near-ground rotation is still generated at 0259 UTC. At 200 m AGL, this rotation has a vertical vorticity greater than 0.1 s−1, reaching the mature stage from 0302 to 0308 UTC with a vertical vorticity greater than 0.2 s−1. The rotation dissipates at 0315 UTC when it approaches the mountainous area.
In WRFLES50m, a strong rotation (tornado) is successfully generated at 0253 UTC, which is 6 min earlier than the generation of the near-ground rotation in WRFLES250m. However, the two rotations are still regarded as physically consistent because the near-ground vorticity development processes are morphologically similar. Compared with the actual tornado, the simulated tornadogenesis occurs about half an hour earlier and the simulated path deviates by approximately 12 km north (Fig. 1b). Other studies of this tornado also show earlier tornadogenesis and north deviation of the path (Mashiko 2016a; Seko et al. 2015). The hydrometeor distribution of the supercell (Fig. 3a) shows a similar pattern to the observations (Yamauchi et al. 2013) and M16 (see his Fig. 2). A low pressure region is located within the hook echo of the supercell and extends to the tornado center. The central pressure deficit of the tornado near the ground is approximately 25 hPa (Fig. 3b).
The time series of the maximum vertical vorticity and maximum horizontal wind speed at the lowest scalar model level (approximately 5 m AGL; ζz=5m, VHz=5m, respectively) are plotted in Fig. 4. ζz=5m and VHz=5m are relatively small during the first 20 min of WRFLES50m. ζz=5m increases quickly from 0252 UTC and exceeds 1.2 s−1 at 0254:30 UTC, while VHz=5m increases until 0255 UTC and exceeds 60 m s−1. The period between 0253 and 0255 UTC is indicated as the tornadogenesis and intensification stages. The minimum sea level pressure drops by 25 hPa over this period. From 0255 to 0310 UTC, ζz=5m and VHz=5m fluctuates between 0.45–1.1 s−1 and 43–63 m s−1, respectively. Starting at approximately 0310 UTC, the decay stage can be indicated by the decreasing tendencies of ζz=5m and VHz=5m.
4. Methodologies: Circulation and vorticity analyses
To quantitatively investigate the vorticity sources, circulation and vorticity analyses are applied on material circuits (attached to air parcels) and air parcels from the near-ground vortices, respectively. The locations of the air parcels are calculated by backward trajectories.
The backward trajectories are calculated using velocity data of WRFLES50m saved every 0.5 s (two model time steps). A fourth-order Runge–Kutta (RK4) method is used with a time step of 0.5 s. The velocity at the air parcel location is obtained by tricubic and linear interpolation in space and time, respectively. Dahl et al. (2012) pointed out that even when the model output saved every model time step are used, trajectories that are traced backward from tornadoes could still be spurious. In the present study, both the model and backward trajectory utilize smaller time steps than most recent studies of the same model grid size (Schenkman et al. 2014; M16; Y18). To check the sensitivity to the saving frequency of model output, the circulation analysis of a tornadic vortex is also tested based on model output saved every four model time steps (not shown). In most periods of the analysis, the difference in the circulation budget is slight between the case of every four model time steps and the default (every two model time steps) case. Nonetheless, the authors still note the caveat that the backward trajectory calculations may result in spurious trajectories.
Brioude et al. (2012) mentioned the uncertainties in interpolating the wind velocity field from the Eulerian model to the Lagrangian model. The authors suggested avoiding using the Cartesian vertical velocity when the terrain is complex. In this study, the vertical velocity is represented by the instantaneous mass-weighted sigma dot vertical velocity (). The is defined on the vertical coordinate of the WRF Model and is directly used in the advection terms of the governing equations (Skamarock et al. 2008). Thus, using is less error prone than using the Cartesian vertical velocity. Meanwhile, the terrain of the regions related to backward trajectory analyses is quite flat (see the Kanto Plain in Fig. 1). Thus, the analyses in this study should not be subject to large error induced by complex terrain.
When the air parcel is below the lowest scalar model level, there is uncertainty in specifying the location, velocity and scalar values of the parcel. Guchte and Dahl (2018) explained the unphysical tendency of the scalar variables on the trajectory below the lowest scalar model level in a simulation using free-slip bottom conditions. For the semislip conditions, as in the present study, the flow is still poorly defined below the lowest scalar model level, and it may lead to errors in the trajectory analyses. Previous semislip studies assumed a logarithmic wind profile (Mashiko et al. 2009; M16), while in the present study, the horizontal velocity is determined by the second-order polynomial extrapolation from the upper levels. Using the extrapolation, the interpolated vorticity and the integrated vorticity (see section 6) show better consistency than those calculated with the logarithmic profile. Y18 also mentioned that using the linear extrapolation results in more consistent circulations than using the logarithmic profile. The authors suggested that the stratification in the tornado was not neutral, and the logarithmic profile might not be suitable for tornadic conditions.
Here ωs, ωc, ζ represent the three components of vorticity under the seminatural coordinates (s, c, z),where s, c, z are horizontal streamwise, horizontal crosswise, and vertical coordinates, respectively. In the present study, the streamwise (and crosswise) vorticity is defined with respect to the ground-relative wind unless otherwise specified. The vorticity equations of an air parcel are given as (Mashiko et al. 2009; Schenkman et al. 2014; Roberts et al. 2016; Y18)
where ψ = tan−1(υ/u) represents the direction of the horizontal ground-relative wind VH = (u, υ); w, p, ρ represent vertical velocity, pressure, and density, respectively; and Fs, Fc, Fz represent three components of the diffusion terms, which consist of subgrid-scale (SGS) mixing and implicit terms. The SGS mixing terms are the spatial derivative of momentum stress. The momentum stress consists of the surface shear stress (parameterized in the surface layer scheme) and the internal mixing stress (the momentum stress above the lowest level, calculated by the LES model). The implicit terms are from the upwind advection schemes. At the lowest two levels, the vertical advection uses central difference; thus, the implicit terms are zero. For the air parcels analyzed in the present study, the implicit term is much smaller than the SGS mixing term. This result is likely because most air parcels are located at low levels. The first three terms on the right-hand side (rhs) of the equations are the stretching and tilting terms, which embody the effects of vortex deformation with fluid. The fourth and fifth terms on the rhs are the baroclinic (solenoidal) and frictional terms. The last terms on the rhs of Eqs. (1) and (2) represent the conversion between crosswise vorticity and streamwise vorticity as the parcel changes moving direction.
In calculating the terms of the vorticity equations, errors may be induced by the large velocity gradient in a highly converged flow (Dahl et al. 2012). Circulation analysis can be used to alleviate this problem by distributing the material circuit surrounding the vortex. The circulation C is defined as
in which υ = (u, υ, w), and l is the line vector along the circuit.
The time change of the circulation can be derived as
The circuit is traced back by the backward trajectories of many air parcels along the circuit. To exclude ambiguity, the initial times of the circuits mentioned in this study refer to the times when the backward trajectory calculations start. Initially, air parcels are distributed along the circuit at 1-m intervals. During the backward integration, if the distance between two adjacent air parcels exceeds 50 m, an additional parcel is added between the two.
The source terms in the circulation and vorticity equations [rhs of Eqs. (1)–(3) and (5)] are interpolated from the nearby model grid points to air parcel locations. If an air parcel falls below the lowest scalar level, it is uncertain to value the source terms, which may lead to errors in the analysis. In the present study, the source terms at the lowest scalar level are used (i.e., constant value). Calculation of the source terms using second-order extrapolation from the upper levels has also been tested. The sum of the time-integrated source terms on the rhs of Eq. (5) is more consistent with the C(t) in Eq. (4) if using the constant value (not shown).
5. Circulation analysis
a. Tornadogenesis process
Figure 5 shows the evolution of the vertical vorticity, horizontal velocity, and vertical velocity at 200 m AGL from 0250 UTC to tornadogenesis (0253 UTC). The FFGF and RFGF are indicated by the regions of horizontal velocity convergence. The FFGF is regarded as the left-flank convergence boundary (LFCB) described by Beck and Weiss (2013) because its orientation and location are different from those of a traditional FFGF (Lemon and Doswell 1979). The updraft region is aligned with the FFGF and RFGF, which is especially strong near the region where the FFGF intersects with the RFGF (hereafter, intersection zone). At 0250 UTC, two vortices exist (Fig. 5a-1). One vortex is located near the intersection zone. This vortex is regarded as the pretornadic vortex owing to the consistency with the tornado (Figs. 5a-2–5a-4). Another vortex is approximately 2 km to the southwest of the pretornadic vortex. This vortex is connected with a primary low-level mesocyclone (LMC, defined below the 1-km height);2 thus, it is named the primary near-ground mesocyclone (NMC). Several local vertical vorticity extrema exist in the primary NMC. The vorticity extrema also evolve, rotating around each other (see online supplement). The RFD region is located behind (on the south and west sides of) the RFGF. In the RFD region, flow of locally intensified horizontal momentum is evident (Fig. 5b-1). This flow is from a locally intensified downdraft, which corresponds to an RFD surge. [The RFD surge reaches the ground around 0250 UTC; thus, the downdraft is less evident in the panels of vertical velocity (Fig. 5c-1)]. Horizontal convergence is enhanced at the front edge of the RFD surge. At 0250 UTC, the RFD surge is located on the south side of the primary LMC and reaches the ground to the south of the intersection zone. Then, the RFD surge heads quickly northeastward, approaching the intersection zone (0251 UTC; Figs. 5b-2,c-2). When the RFD surge reaches the RFGF (0252 UTC; Figs. 5b-3,c-3), the updraft on the left front edge of the RFD surge is merged with the updraft in the pretornadic vortex, leading to the intensification of vertical vorticity (tornadogenesis at 0253 UTC; Figs. 5a-4,b-4,c-4). This indicates that the RFD surge may play a role in triggering tornadogenesis by increasing the convergence near the ground. The phenomenon that an RFD surge causes a tornado to the east of a primary mesocyclone is similar to cyclic mesocyclogenesis (Adlerman et al. 1999; Adlerman 2003; Beck et al. 2006).
To be more intuitive, the evolution of the primary mesocyclone and tornadogenesis is shown by 3D rendering in Fig. 6. Before tornadogenesis (Fig. 6a) the primary LMC and NMC are connected to an MMC (named the primary MMC). A new LMC is formed at approximately 0251 UTC (Fig. 6b) above the pretornadic vortex. During the intensification of the new LMC, the pretornadic vortex is rapidly intensified and connected to the new LMC (Fig. 6c). Eventually, the pretornadic vortex becomes the tornado after 0252 UTC (Fig. 6d). The new LMC is named the tornadic LMC, hereafter.
To illustrate the changes in the flow entering the near-ground vortex at approximately the onset of tornadogenesis, circuits initialized from the pretornadic and tornadic vortices are traced backward for 10 min (Fig. 7). The spatial coverages of the circuits after backward integration could qualitatively represent the distributions of air parcels that were initially enclosed by the circuits. The circuits are initiated as rings3 surrounding the cores of the vortices at 200 m AGL (e.g., Fig. 8a). Four initial times are selected from 0250 to 0253 UTC. For the circuit initiated at 0250 UTC, most part of the circuit is traced back near the surface in the eastern inflow.4 For the circuits initiated after 0251 UTC, segment of the circuits are traced back to the FFGF region. The circuit segment in the FFGF region is gradually traced back to approximately the 300-m height, which means the segment has experienced a downdraft. Considering that the FFGF lies between the cold pool and the warm updraft, the horizontal vorticity could be generated by baroclinity. When an air parcel traveling along the FFGF experiences a downdraft, it leads to a “slippage” between the velocity and vorticity vectors. Thus, the horizontal vorticity is tilted into vertical vorticity during the descent, and then, the vertical vorticity can be intensified by stretching (Davies-Jones and Brooks 1993; Markowski and Richardson 2014; Rotunno et al. 2017). For the circuits initialized after 0252 UTC, large segments of the circuits are traced back to the RFD region on the south and west sides of the primary LMC/NMC. It is reiterated that the RFD surge may play a role in triggering tornadogenesis. Note that the segments of circuits become very contorted. The rotational flow and the local vorticity extrema in the primary LMC/NMC may cause this contortion (see online supplement of the detailed evolution in a zoom-in view). As suggested by Dahl et al. (2012), a strongly curved, confluent flow may induce errors in the backward trajectory calculations because the gradients of the velocity field could be poorly resolved. A caveat is noted here about the possible errors in backtracking the circuits.
To investigate the change in the circulation sources before and just after tornadogenesis, the circulation sources of the near-ground vortices at the pretornadic (0250 UTC) and tornadogenesis (0253 UTC) times are discussed in the following subsections.
b. Circulation analysis of the pretornadic vortex
The circuit initialized surrounding the pretornadic vortex is shown in Fig. 8. The circulation calculated by Eq. (4) is compared with the circulation integrated from the terms on the rhs of Eq. (5) (named the integrated circulation). The histories of the two circulations show good agreement, which verifies the reliability of the circulation calculation. The magnitude of the frictional term is much larger than that of the baroclinic term. The positive frictional term lasts from approximately 0230 to 0234 UTC. After 0234 UTC, the frictional term is negative, resulting in a decrease in circulation. The decrease in circulation just before a circuit arrives at a vortex is also noticed in previous observational (e.g., Markowski et al. 2012) and numerical (e.g., M16; Schenkman et al. 2012) studies. The circuit is located near the ground (mostly lower than 200 m AGL; see Fig. 7a), implying that the frictional force reduces circulation when the flow converges near the ground.
c. Circulation analysis of the vortex at tornadogenesis time
The circulation analysis of the vortex at tornadogenesis time is shown in Fig. 9. The frictional term is also the dominant source. The frictional term is positive throughout most of the period. The negative frictional term also occurs during the last 1.5 min before the circuit arrives at the tornado. The baroclinic term is negligible compared to the frictional term. This result apparently differs from M16, who emphasized the contribution of baroclinity.
According to the caveat mentioned in section 5a, the complex flow in the primary LMC/NMC may induce error in the trajectory analysis. In addition, it could be supposed that if the flow is more complex, the error in trajectory calculation are more sensitive to the saving frequencies of model output. Evidence for the error is shown in the sensitivity tests. Prior to min 15 in Fig. 9, the period in which large circuit segments are contorted by the primary LMC/NMC, discernable differences in the circulation and the frictional term appear between the default case (model output saved every two time steps) and the case of every four time steps (not shown). The errors in trajectory calculation may also be the reason for the intensive oscillation of the frictional term prior to min 15.
The configurations of the circuit traced back in time are shown in Figs. 10a–10e. The aforementioned circuit segment on the south and west sides of the primary LMC is present before 0249 UTC. A primary hook echo is associated with the primary LMC. M16 found that the hook echo was accompanied by large baroclinity, which was the vorticity source of the tornado. To investigate the weak baroclinity in the present study, the buoyancy around the primary hook echo at 0249 UTC is plotted (Fig. 10f). The buoyancy term is estimated from the Boussinesq version of the vertical momentum equation:
where g is the gravitational acceleration, ρ′ is the density deviation from , p′ is the perturbation pressure, which is the deviation from the base-state pressure, and is the average density over the 7.5 km by 7.5 km horizontal region around the hook echo. The first and second terms on the rhs are the vertical perturbation pressure gradient force (VPPGF) and buoyancy terms, respectively.
The buoyancy term difference across the hook echo is approximately 0.06 m s−2, which is smaller than that in M16 (approximately 0.2 m s−2, see his Fig. 11b). This difference between the two studies implies that compared with M16, the present study has weaker baroclinity at the hook echo and a warmer cold pool. The microphysics schemes may cause the difference. M16 utilized a single-moment scheme in which a fixed intercept parameter (8 × 106 m−4) was used to estimate the raindrop size distribution. This value is typically used in a number of previous studies (e.g., Schenkman et al. 2014), but there is still uncertainty regarding this parameter which is adjusted in some other studies (e.g., Dahl et al. 2014; Roberts et al. 2016). However, under realistic conditions, the intercept parameter of raindrops varies over a wide range. Multimoment schemes diagnose both mixing ratios and number concentrations; thus, the intercept number does not need to be predefined. Dawson et al. (2015, 2016) also suggested that multimoment microphysics schemes show advantages over single-moment schemes in predicting tornado track and intensity, because single-moment schemes tend to produce cold pools that are too strong. The two-moment microphysics scheme used in the present study supposedly resolves a more realistic condition.
Moreover, because the circuit segment is on the south and west sides of the hook echo rather than across the tip of the hook echo, there is no baroclinic vortex line penetrating the circuit, resulting in the small baroclinic term in circulation budget.5 This configuration of the circuit segment near the hook echo may result from the relationship between the primary LMC/hook echo and the tornado. The primary LMC/hook echo provides the downdraft6 (the RFD surge), which enters the tornado at genesis time. However, it is the tornadic LMC (associated with a tornadic hook echo) that provides the updraft for tornadogenesis. In M16, the circuit connects the downdraft and updraft regions associated with the same hook echo, allowing the baroclinity to be maintained. In the present study, the baroclinity is less likely to be maintained.
The discussion above indicates that the baroclinity at the hook echo is not necessarily a major vorticity source responsible for tornadogenesis. Here, we reiterated that frictional term is the major vorticity source. The fractions of the air parcels along the circuit in different height bins are shown in Fig. 11a. After being traced back for 3 min, most of the air parcels are below 100 m AGL. The fraction of the air parcels near the ground increases when the circuit is traced backward to an earlier time. Physically, this implies that a large proportion of the air in the tornado is drawn from the near-ground levels. The distribution of frictional term (Fig. 11b) is similar to that of the fraction of air parcels. The most contribution from friction occurs near the ground, especially in the first layer (below 10 m AGL). This result implies that surface drag has the dominant contribution to the frictionally generated circulation.
6. Vorticity budget analyses of air parcels
To explore the large frictionally generated vorticity in detail, the vorticity budgets of 140 air parcels in the tornado at genesis time (0253 UTC) are analyzed (Fig. 12a). The air parcels are initialized at 200 m AGL and are located at the centers of the 140 grids enclosed by the circuit analyzed in section 5c. The trajectories of the air parcels are calculated backward for 10 min (Fig. 12b). Based on where the air parcels originate, the air parcels can be categorized into three classes: (class-A) air parcels descending with the RFD surge, (class-B) air parcels from the north, and (class-C) air parcels from the low level in the eastern inflow.7
All three classes of the air parcels have contributions to the total circulation (the area integral of the vertical vorticity) of the tornado (Fig. 12c), while class-C contributes the most (60.6%). Class-C also occupies 44.3% of the total number of air parcels (Fig. 12d). Class-A and class-B only contribute 10.4% and 29% of the total circulation, respectively. This finding indicates that most vorticity sources should be attributed to the eastern inflow (class-C) and may also indicate that the main function of the RFD surge (class-A) is to enhance the near-ground convergence rather than to provide vorticity. It should be noted that these indications are only based on the tornado at genesis time. The conclusions could be different in other periods of the tornado (Marquis et al. 2016).
For trajectories in each class, histories of the vorticity terms show qualitatively similar tendencies. The vorticity calculated by interpolation from the vorticity on nearby model grid points is compared with the vorticity calculated by integrating the terms on the rhs of the vorticity equations. The former vorticity is named the interpolated vorticity. The latter is named the integrated vorticity. All of the trajectories show good agreement between the two vorticities (look ahead to Figs. 13c, 14c, and 16c for examples). This result indicates that the calculated trajectories and vorticity terms are reliable. In each class, the air parcel that has the best agreement is selected as a representative to be discussed in detail.
a. Class-A: Air parcel descending with the RFD surge
The representative air parcel of class-A (Fig. 13) starts descending from approximately 200 m AGL at 4 min before tornadogenesis (T − 4 min, hereafter). From T − 2 min to approximately T − 0.25 min, the air parcel is near the ground and moves quickly northeastward toward the tornado (also described in section 5a). During this period, crosswise vorticity increases mainly via frictional term. Streamwise vorticity increases mainly via the exchange from the crosswise vorticity. As mentioned previously, the streamwise and crosswise directions are defined with respect to the ground-relative wind. Just before arriving at the tornado (from T − 0.25 min), the air parcel moves northward in the east of the updraft region. Therefore, the crosswise and streamwise vorticity are tilted into the vertical (Fig. 13a). After the air parcel enters the tornado, the increase in vertical vorticity is dominated by stretching.
b. Class-B: Air parcel from the north
The air parcels in class-B can be further divided into two groups. In the first group the air parcels originate from the cold pool to the west of the FFGF (green color in Fig. 12), while in the second group the air parcels are along the FFGF (purple color in Fig. 12). A distinct feature of the first group is that the air parcels experience several “internal” boundaries to the west of the FFGF. The representative air parcel is from the first group (Fig. 14). Before arriving at the FFGF, the air parcel transverses three “internal” boundaries (Fig. 15; as mentioned in section 5a, the FFGF shows consistency with the LFCB, thus labeled as “LFCB”). Updrafts are associated with the “internal” boundaries and there are downdrafts between the boundaries, which cause the air parcel to oscillate vertically.
From T − 0.5 to T − 0 min, the tilting of crosswise vorticity is predominant for the increase in the vertical vorticity. The crosswise vorticity is generated mainly owing to frictional and tilting terms from T − 1.5 to T − 1 min when the air parcel is at its nadir. This observation indicates that the frictional term may result from surface drag. Schenkman et al. (2014) generally found that large frictional generation was present when the parcels were closer to 50 m AGL and lower, while herein, the nadir of the air parcel is approximately 60 m AGL thus the frictional term is relatively small. For several other air parcels in the first group that have nadirs lower than 60 m AGL, the frictional terms are more significant (not shown).
To investigate the vertical influential range of the surface drag, the history of the vertical profile of the crosswise frictional term8 at the location of the representative air parcel is checked (not shown). It is found that from T − 1.5 to T − 1 min, the frictional term is enhanced between 50 and 150 m AGL. This enhancement is probably caused by the strong updraft in the FFGF where the air parcel is located (Figs. 15a-4,b-4). The strong updraft may extend the vertical range of the frictional term. At the levels above the lowest one, more active internal mixing stress may be induced by the large velocity gradient associated with the updraft and the gust front, and consequently results in the large frictional term.9 Similar (but weaker) enhancement is also found on the “internal” boundaries, which is consistent with the (weaker) updraft and convergence on these boundaries.
For all the air parcels in class-B, the baroclinic term is quite small, which may result from the weak density gradient across the FFGF. The potential virtual temperature difference between the cold pool and the warm inflow is approximately 4 K in this study (not shown). In M16, this difference is much larger, which is approximately 15 K (W. Mashiko 2018, personal communication). The result of the present study is different from those in several previous studies that emphasized the baroclinic effects (Wicker and Wilhelmson 1995; Markowski and Richardson 2014; Dahl 2015). In the second group, the air parcels originate from lower levels. The large crosswise vorticity is generated by friction. The baroclinic terms are also insignificant. The histories of the air parcels are similar to those of the air parcels in class-C except for the different moving directions; thus, a detailed description is omitted.
c. Class-C: Air parcel from the low level in the eastern inflow
The ground-relative inflow of the “eastern” parcels is actually from the south. Thus, the eastern inflow is also named southerly flow in this subsection. From T − 9 to T − 2 min, the crosswise vorticity of the representative air parcel (Fig. 16) increases owing to frictional term. The air parcel is near the ground during this period. Figure 17a depicts the history of the vertical profile of the crosswise frictional term at the location of the air parcel, showing that the discernible magnitude of the term is mainly below 50 m AGL. This implies that the frictional term is mainly induced by surface drag. From T − 0.5 to T − 0 min, the air parcel enters the tornado from the east side of the updraft. Thus, the crosswise vorticity is tilted into the vertical because the updraft region is located on the left side of the air parcel motion.
The environmental crosswise vorticity is also quite large (see T − 10 min in Fig. 16c). Because the environmental inflow (southerly inflow) is almost steady (see the eastern part in Fig. 16a-1), it is thought that the environmental inflow is similar to the flow in a typical boundary layer. Additionally, there is no comparable crosswise vorticity above 100 m AGL in the environmental inflow (see T − 10 min in Fig. 17b). Thus, this result indicates that the environmental crosswise vorticity (vertical wind shear) is also caused by surface drag.
Although the vertical grid is fine near the ground (10-m spacing), several air parcels in class-C are still traced back below the lowest scalar model level. Higher vertical resolution near the ground is necessary in future research (Rotunno et al. 2017).
7. Discussion and conclusions
As mentioned previously, for the representative air parcels, the tilting of the crosswise vorticity into the vertical is predominant. However, significant upward tilting of streamwise vorticity is also found in several other parcels in each class. The streamwise vorticity is attained via exchange from crosswise vorticity when the moving directions of the air parcels change. This behavior is consistent with the second mechanism demonstrated by Roberts et al. (2016), and it is more prominent for the air parcels that have experienced large changes in moving directions (e.g., the black-marked air parcels in the south of the vortex in Fig. 12a, which was also noted in Markowski 2016). Furthermore, after the representative air parcels enter the tornado and swirl, the streamwise vorticity would increase and its tilting into vertical vorticity would be more significant.
Davies-Jones (1984) demonstrated that to form the updraft rotation (helical vortex) in supercell storms, the environmental horizontal vorticity should have a streamwise (in storm-relative frame) component (see his Figs. 7,8). This concern may raise some doubt about whether the tilting of crosswise vorticity in the present study can truly turn into the helical vortex of the tornado. Actually, as mentioned previously, the “crosswise” in the present study is ground-relative.10 We checked the vorticity in storm-relative frame and found that the tilting of storm-relative streamwise vorticity is prominent (Fig. 18). Figures 18a–18d show the ground-relative and storm-relative flow and the vortex lines passing the representative air parcels during (T − 10 s) the tilting of horizontal vorticity. The horizontal vorticity vectors are roughly parallel (perpendicular) to storm-relative (ground-relative) horizontal velocity vectors (especially evident in the eastern inflow region). The storm-relative streamwise vorticity is turned into the vertical (Figs. 18b,d–g) and the helical vortex is formed. Thus, the tilting of horizontal vorticity in this study is consistent with Davies-Jones (1984).
In this study, the 6 May 2012 Tsukuba tornado was successfully reproduced utilizing WRFLES. The tornado was explicitly resolved with 50-m horizontal grid spacing. The vorticity sources responsible for the tornadogenesis were investigated by circulation and vorticity analyses. The following conclusions can be drawn.
Circulation budgets for the circuits initialized around the pretornadic and tornadic vortices showed that the magnitude of the frictional term was consistently much larger than the magnitude of the baroclinic term. In the case of the pretornadic vortex several minutes before tornadogenesis, friction actually reduced the circulation when the circuit approached the vortex. For the vortex at the time of tornadogenesis, friction contributed significant positive circulation. Considering the controversy of the vorticity sources of tornadoes, the present study emphasized again the dominant contribution of frictionally generated vorticity as proposed by Schenkman et al. (2014). The air parcels entering the tornado at genesis time could be categorized into three classes: (class-A) air parcels descending with the RFD surge, (class-B) air parcels from the north, and (class-C) air parcels from the low level in the eastern inflow. For each class, the surface drag was the main vorticity source of the tornado. The frictionally generated crosswise vorticity was directly tilted (or first exchanged into the streamwise vorticity and then tilted) into the vertical. Among the three classes, class-C contributed the most to the circulation of the tornado.
An RFD surge from the downdraft region on the south and west sides of the primary LMC/hook echo might play a role in triggering a tornado by increasing the convergence near the ground. The baroclinity across the hook echo was weaker than that in M16, which might lead to the lack of baroclinity in the vorticity sources.
Caveats should be mentioned that the MOST surface layer scheme might be questionable in the outflow; the LES might suffer from the overestimated wind shear due to underresolved turbulence and inadequate grid resolution near the surface. Another caveat was that as the air parcels were allowed to descend below the lowest scaler model level, the backward trajectories were prone to errors. The errors consisted of the uncertainty in both the location of air parcels and the extrapolated relevant fields. Moreover, the backward trajectory calculations might still need a longer time to reveal all the possible vorticity sources, especially considering that after being calculated backward for 20 min, the circulation of the circuit from the pretornadic vortex was even larger than the initial circulation. The circulation of the circuit from the tornado almost reached zero after a 13-min backward calculation. However, the circulation source terms in a prior time might still be uncertain. This uncertainty can be indicated from a circuit studied in Markowski (2016). The baroclinic term was not prominent until the circuit was calculated backward for 15 min, and the frictional term showed a remarkable negative value after a 25-min backward calculation (see his Fig. 23).
Dr. Hidenori Kawai contributed a considerable amount of help in the simulations of this study. The authors appreciate the valuable advice about the numerical configurations and trajectory calculations from Dr. Wataru Mashiko. The authors also acknowledge Prof. Hiroshi Niino, and Drs. Hiromu Seko and Kazuo Saito for the helpful discussions about the WRF simulation of tornadoes. We also thank Christopher Weiss and the anonymous reviewers for their constructive comments and suggestions on the original manuscript. The simulations were performed using the supercomputer TSUBAME 3.0 at the Tokyo Institute of Technology, financially assisted by the MEXT FLAGSHIP project: “Advancement of meteorological and global environment predictions utilizing observational ‘Big Data’” of the social and scientific priority issues (Theme 4). The authors also appreciate the support from the China Scholarship Council (CSC) for overseas doctoral study.
The horizontal grid spacing in their simulation was 500 m. Thus, the tornado was not explicitly resolved. The near-surface rotation that presented characteristics of a tornado is referred to as near-surface mesocyclone.
In WRFLES250m, the LMC is defined where vertical vorticity (ζ) at the 1-km height is larger than 0.05 s−1. ζ is concentrated at the intersection between FFGF and RFGF; thus, the LMC is easy to identify. In WRFLES50m, at the corresponding location there are several local vorticity extrema.
Circuits located along vorticity isopleths have also been tested. The resulting circulation budgets and locations are qualitatively similar to those of circuits initiated as rings.
In the storm-relative frame, the inflow direction is westward and is thus named the eastern inflow.
The downdraft is mainly forced by the negative VPPGF term on the rear side of the primary LMC (not shown).
About five air parcels in class-C (trajectories from the southmost in Fig. 12b) are traced back to a relatively higher location in weak convection from the south. These air parcels also have a similar history as the other air parcels in class-C except for some pre-existing vertical vorticity.
In Dahl et al. (2012) the air parcels that are traced back to the eastern inflow are regarded as erroneous parcels because the storm outflow forms a boundary preventing the inflow from entering the tornado near the ground. The authors utilized free-slip bottom condition. On the other hand, several recent studies did show air parcels coming from the eastern inflow (e.g., Naylor and Gilmore 2014; Schenkman et al. 2014; Marquis et al. 2016; M16; Roberts and Xue 2017). Most of these studies used nonfree-slip bottom condition. This may indicate that the nonfree-slip condition alters storm outflow features, enabling the eastern inflow to enter the low level of the tornado. Despite this, the possible error in trajectory is noted here.
This explanation may lead to a thought that internal mixing is also a vorticity source. However, Markowski (2016) found that in free-slip simulation, internal mixing actually decreases the circulation. Thus, we consider that internal mixing should not be regarded as a vorticity source. Surface drag is the vorticity source, and the effect of surface drag is extended upward via internal mixing.