Investigation into the Formation, Structure, and Evolution of an EF4 Tornado in East China Using a High-Resolution Numerical Simulation

Devastating tornadoes in China have received growing attention in recent years, but little is known about their formation, structure, and evolution on the tornadic scale. Most of these tornadoes develop within the East Asian monsoon regime, in an environment quite different from tornadoes in the U.S. In this study, we used an idealized, highresolution (25-m grid spacing) numerical simulation to investigate the deadly EF4 (Enhanced Fujita scale category 4) tornado that occurred on 23 June 2016 and claimed 99 lives in Yancheng, Jiangsu Province. A tornadic supercell developed in the simulation that had striking similarities to radar observations. The violent tornado in Funing County was reproduced, exceeding EF4 (74 m s–1), consistent with the on-site damage survey. It was accompanied by a funnel cloud that extended to the surface, and exhibited a double-helix vorticity structure. The signal of tornado genesis was found first at the cloud base in the pressure perturbation field, and then developed both upward and downward in terms of maximum vertical velocity overlapping with the intense vertical vorticity centers. The tornado’s demise was found to accompany strong downdrafts overlapping with the intense vorticity centers. One of the interesting findings of this work is that a violent surface vortex was able to be generated and maintained, even though the simulation employed a free-slip lower boundary condition. The success of this simulation, despite using an idealized numerical approach, provides a means to investigate more historical tornadoes in China.


Introduction
Tornados are a devastating hazard that occurs at a very small scale and has a short lifetime. Direct observations of tornadoes are difficult to obtain (Wurman et al., 2013) and require in-situ campaigns (Wurman et al., 2012). Because of this, numerical simulation has emerged as an approach to investigating  and predicting tornadoes . Although progress has been made in numerical simulations of severe weather, the reproduction and prediction of small-scale tornadoes remain challenging (Markowski and Richardson, 2009). One reason is the lack of calculation resources, as the limited resolution of numerical models cannot truly resolve the physical processes operating close to the scale of a tornado (Sun et al., 2014). The other is our poor understanding of, or inability to reproduce, these crucial physical processes (Stensrud et al., 2013); hence, the processes resolved by high-resolution models may be highly artificial or incomplete. For instance, even the initiation of convection has not yet been satisfactorily resolved, let alone the genesis of a tornado.
As an alternative, a semi-idealized simulation approach could be employed, in which a high-resolution cloud model and a base-state sounding profile are combined (Coffer and Parker, 2017). This approach is referred to as "semi-idealized" because the initial conditions are horizontally homogeneous and initialized by artificial initiation, but the sounding used to generate the initial conditions is obtained from real case data-quite different from most previous studies, which have tended to be based on idealized sounding profiles also (Davies-Jones and Markowski, 2013). The base-state sounding can be taken from near-storm observations or real-case numerical weather prediction (NWP) models, and an artificial initialization of thunderstorm convection can be used in the cloud model. This approach was first used to investigate supercell formation, enhancement, and maintenance, at the relatively low resolution of around 1 km, as well as warm rain microphysics (Klemp et al., 1981;Grasso and Cotton, 1995;Wicker and Wilhelmson, 1995). In the past few decades, with the growth of computing resources, numerical simulations with resolutions of around 100 m have been performed, with more sophisticated double-moment microphysical schemes; however, such simulation studies have focused mostly on the circulation and origins of strong surface vortices, rather than a real tornado characterized by a funnel cloud, since the grid spacing was still too coarse Naylor and Gilmore, 2014;Nowotarski et al., 2014;Markowski and Richardson, 2017;Rotunno et al., 2017). Numerically simulated tornadic storms that adequately reflect their observed counterparts are quite rare in the literature. It was not until very recently that the detailed structure of a simulated tornado (with a grid spacing of 30 m) was formally reported (Orf et al., 2017). Considering the complexity and current knowledge limitations involved, the numerical simulation of tornadoes is a hot topic that awaits further study worldwide.
In China, tornadoes have not been a particularly wellstudied weather hazard, partly because of their low frequency of occurrence compared with other weather systems that strike more frequently and with greater impacts, such as typhoons, heavy rainfall, flooding, and drought. In truth, they are quite common (Yao et al., 2015), and have thus started to cause more concern in recent years. Damage surveys of high-impact tornado events have been performed in recent years to obtain more detailed information about tornado or possibletornado events, such as the dangerous winds during the heavy rainfall event that occurred in Beijing on 21 July 2012 (Meng and Yao, 2014); the overturning of the Oriental Star in Jianli, with over 450 fatalities (Meng et al., 2016); the tornadic supercell during the landfall of Typhoon Mujigae in 2015 (Bai et al., 2017;Zhao et al., 2017); and several others . These studies put forward a series of developments in our understanding of the structure and formation of tornadoes in China, but the ability of damage surveys and radar observations alone to provide a complete life story of tornadoes is limited. Numerical simulation of the structures and processes on the scale of tornadoes is urgently needed.
In the present study, we began by assessing the feasibility of a method to simulate tornado events. Then, after successfully reproducing an EF4 (Enhanced Fujita scale category 4) tornado that occurred in Funing, Yancheng, Jiangsu Province, China, we used the simulation to investigate the tornado's formation and evolution, as well as some preliminary characteristic structures and processes that merit further analysis.

Methodology and model configuration
The numerical simulation of the tornado was performed by using Cloud Model version 1 (CM1), release 18.3 (Bryan and Fritsch, 2002). CM1 is a nonhydrostatic cloud model developed and maintained by Dr. George Bryan from the NCAR. In general, for simulations of tornadic supercells, the initial condition used in CM1 is a homogenous background based on a sounding profile. This profile can be either an idealized one that simplifies the vertical distribution of parameters, or an observed one that represents characteristics of the actual conditions. The convection in CM1 is initiated with artificial perturbation, e.g., an ellipsoid warm bubble, forced convergence, etc.
Our simulation used a stretched grid domain to represent finer-scale structures, using a resolution as high as 25 m within the focused area, and to save on computational resource. The vertical resolution was set evenly to 20 m below 1 km AGL, which set the lowest model level to 10 m AGL. The model top height was 18 km. The updraft nudging technique (Naylor et al., 2012) was utilized for the convection initiation. Other configurations of the model were set following previous studies with similar aims. The large time step was set to be 1 s, and a selfadaptive option was chosen. The simulation was run for 3 h, with the output saved every 10 s. A fifth-order advection approximation was used both horizontally and vertically. The pressure solver used was Klemp-Wilhelmson time-splitting for acoustic modes, with a vertically implicit solver, and horizontally explicit calculations, as in the Weather Research and Forecasting (WRF) Model and the Advanced Regional Prediction System (ARPS). The turbulence closure used was the Turbulent Kinetic Energy (TKE) scheme. No topography was included. The explicit moisture scheme used was the Morrison double-moment microphysics scheme, which is more sophisticated and generally produces better simulations than warm-rain schemes. The lateral boundary conditions were set to be open-radiative, while both the bot-tom and the top boundary conditions were free-slip. The free-slip boundary conditions were those of Klemp and Wilhelmson (1978), in which the shear stress at the boundaries is set to the stress at the first interior grid point, as opposed to zero-stress boundary conditions. The reproduction of a tornado exceeding EF4 in this simulation using the free-slip boundary condition is contrary to the previous belief that the "near-ground cyclone collapses into a tornado only if the model includes surface friction" (Davies-Jones, 2015), which is consistent with the findings of Orf et al. (2017). Although there are some indications that surface friction can be a non-trivial vorticity source for tornadoes Roberts et al., 2016), a free-slip condition was used for our simulation given the considerable uncertainty about how to appropriately account for the effects of surface friction in the lower boundary condition, especially when there is no reason to expect the vertical profile of horizontal wind to obey a "log law", such as in the case of unsteadiness and horizontal inhomogeneity (Markowski, 2016). Moreover, the use of an large eddy simulation (LES) turbulence closure in conjunction with surface drag can also be problematic (the surface vertical shear is greatly exaggerated) if the resolved flow is insufficiently turbulent (Markowski and Bryan, 2016). The sensitivity of the formation and characteristics of the simulated tornado to the lower boundary condition will be explored in a subsequent paper.
The base-state sounding used for the initial condition in CM1 was obtained through a reasonable WRF simulation with grid spacing of 3 km. Version 3.8.1 of WRF was used. The simulation contained three domains, with the innermost domain of 3-km grid spacing centered on the observed tornado track. The horizontal grid length of the outermost domain was 27 km, covering most of the related large-scale meteorological systems in East Asia. There were 86 vertical levels, stretching from roughly 20 m above the surface to thicker layers in the upper troposphere, and then thinner layers at the top of the stratosphere. Multiple physical processes important in capturing the synoptic systems and the development of convective storms were included in the simulation. The dynamics was set to the default configuration. A rapid radiative transfer model including the Monte Carlo Independent Column Approximation (McICA) method (Iacono et al., 2008) was used for the shortwave and longwave radiation. The convective parameterization used was the new eta Kain-Fritsch convection scheme (Kain, 2004), and the microphysics scheme adopted was the WRF double-moment six-class method (Lim and Hong, 2010). A Noah land surface model (Ek et al., 2003) was used for the surface processes, along with the Yonsei University Planetary Boundary Layer scheme (Hong and Pan, 1996). The background and boundary conditions used were the six-hourly final analysis (FNL) data provided by the NCEP. The model was initiated at 1200 UTC 22 June 2016, with a diverse range of in-situ surface observations and soundings contained in the domain assimilated by the WRF 3D-Var module. The sounding used to initiate the CM1 simulation was extracted from the 3-km resolution domain near the storm to represent the environment in which the tornado occurred at 0600 UTC 23 June 2016.
The sounding profile had a convective available potential energy (CAPE) of 3955 J kg -1 and strong lowlevel vertical wind shear, which was favorable for tornadogenesis (Fig. 1a). The storm relative helicity (SRH) APRIL 2018 between 0 and 3 km (0 and 1 km) exceeded 200 m 2 s -2 (120 m 2 s -2 ), which is a good indicator for severe tornadoes (Rasmussen and Blanchard, 1998;Rasmussen, 2003). The hodograph plot of the sounding profile showed a veering pattern, which favored tornadogenesis (Fig. 1b). Sensitivity tests showed that, although the general structure of the simulated supercell was robust based on initial soundings, with slight differences, the success or failure of tornadogenesis was highly sensitive to the sounding profile. Among the observations of this case, the closest rawinsonde observation was the 1400 LST sounding at Sheyang, which did not show an environment favorable for tornadogenesis. This strongly suggests that the representativeness of a proximity sounding should be evaluated very carefully. The sounding profile was chosen based on criteria that better favor tornadogenesis. Nevertheless, the most favorable scenarios do not necessarily guarantee representativeness of the real case. Bearing this in mind, our investigation was based on the use of the structures and mechanisms in a plausible proximity environment (and one of the most favorable scenarios) to trigger the tornado, which ensured tornadogenesis and showed some similarities to the observations of the Yancheng EF4 tornado. The sensitivity and selection of a proximity sounding that bears the best resemblance to a particular case await further study and will be presented in a separate paper.

Evolution and structure of the simulated tornado
The CM1 simulation with the selected sounding from WRF successfully captured the tornadic supercell and the embedded tornado. The maximum surface wind in the simulation exceeded 80 m s -1 and lasted for about 10 min, which reached the threshold of an EF4 tornado, which was the rating of this case determined by the damage survey. This section begins by comparing the model results with the observations, and then investigates the evolution and structure of the simulated tornado.

Verification of the simulation by radar observations
The S-band Doppler radar in Yancheng (station number: Z9515, location: 33.4°N, 120.2°E) successfully captured the tornadic supercell (Fig. 2). The general size and morphology of the simulated supercell immediately after tornadogenesis resembled the observations (Figs. 2a,d) at the 0.5° scanning elevation angle, including the hook echo and weak echo region. One intriguing difference between the simulation and the observations was the echo strength of the hook echo. In the simulation, the center of the hook echo had a low value of reflectivity. In contrast, the observations showed high reflectivity echoes. One possible reason for this is the influence of the debris elevated by the strong tornado, which can be observed by the radar as a tornado debris signature (Houser et al., 2016). However, in the current simulation, neither the reflection of the radar signal by the debris, nor the debris itself, was included. Using a tornado simulation with debris algorithms would be a helpful method to solve this mystery (Bodine et al., 2016), but this is beyond the scope of the current study.
The simulated radial velocity at 0.5° PPI (plan position indicator) showed a distinct tornadic vortex signature (TVS) of stronger intensity than the observations. The rotation of the mesocyclone, however, was weaker in the simulation, showing a convergence signature with weaker inbound velocity and stronger outbound velocity towards the Doppler radar (Figs. 2b, e). These features indicate weaker outflow (northwestern wind) from the rear side (west) of the mesocyclone and tornado and stronger inflow (southeastern wind) from the forward side (east). The mid-level mesocyclone was quite similar in both the simulation and observations, as the radial velocity at 4.3° PPI showed mesocyclone structures that were similar in both radius and intensity (Figs. 2c, f). In general, the simulation reproduced the key characteristics of the mesocyclone reasonably, which make it important in revealing the structure and evolution of the tornado, despite differences in terms of the low-level circulations. It should be noted that the differences in the lowlevel velocity field imply that the tornado's structure in the simulation might have been different to that of the actual tornado. Nevertheless, our intention was to carry out a preliminary study on the structure and mechanisms based on a plausible scenario, rather than the unique truth, and so we did not deem this to be a major concern.

Evolution of the simulated tornado
In the 3-h period of the simulation, there were two tornadoes that formed inside the same convective storm. The first tornado was similar to the deadly tornado that struck Funing County and killed 98 people, which was the main focus of the current study. The second tornado may have been representative of a weaker tornado that happened in Sheyang County to the east of Funing during the same event, which caused one fatality. In all, the tornado event on 23 June 2016 killed 99 people in Yancheng City.
The evolution of the simulated storm can be roughly divided into six stages (Fig. 3). It took about 30 min for the updraft-initiated convection to form a storm, and about another 30 min for the storm to evolve into a supercell. After the formation of the supercell, its mid-level mesocyclone intensified and stayed relatively stable (the maximum 2-5 km updraft helicity in Fig. 3a). Multiple surface vortices emerged during this period, but none of them was strong enough to produce a tornado (Fig. 3b). It was not until 90 min into the simulation that the midlevel mesocyclone intensified drastically, followed by the intensification of significant surface vertical vorticity, which indicated the tornado's formation.
The evolution of the tornado was also characterized by negatively correlated growth of the surface velocity and pressure fields (Figs. 3c, d; Fig. 4). At t = 102 min, the maximum surface wind speed grew drastically and reached its maximum at 105 min. Meanwhile, the minimum surface pressure of the tornado dropped from around 980 hPa to below 900 hPa, with the lowest pressure below 850 hPa. The maximum surface wind then weakened for several minutes (though still maintaining the strength of a severe tornado), but soon re-intensified to a second peak at t = 111 min. During this period, the minimum surface pressure fluctuated but generally weakened. This tornado decayed at t = 115 min.
After the demise of the first tornado, although the surface vertical vorticity weakened, the mid-level mesocyclone re-intensified to reach a second peak (Fig. 3a). It then weakened again, but soon reached its third peak at around t = 130 min, followed by re-intensification of the surface rotation, indicating a second tornadogenesis (Fig. 3b). After the demise of this second tornado, the supercell eventually weakened, with no sign of a third tornadogenesis.
The maximum surface vertical vorticity showed two peaks during the lifetime of the simulated storm (Fig. 3b), each of which indicated a concentrated surface vortex for the tornado. Interestingly, the updraft helicity (integrated over 2-5 km AGL) showed three peaks (Fig. 3a). In recent simulation-based studies of tornado events, updraft helicity-defined as the integration of the product of vertical velocity and vertical vorticity at each level over a certain layer (Kain et al., 2008)-has been used to indicate the occurrence of tornadoes and the pos- APRIL 2018 sible tornado tracks (Clark et al., 2012(Clark et al., , 2013. In principle, updraft helicity represents the evolution of rotational updraft intensity in the mid-levels, instead of the near-surface rotation of the tornado. It is more precise to interpret this as the indication of a tornadic mesocyclone.

Structure of the simulated tornado
The supercell in which the tornado was formed had a total height of about 18 km, and it had a cloud base at about 1 km AGL (Fig. 5a). In the forward flank (righthand side of the storm in Fig. 5a, which is the east) of the storm, the cloud base was even and smooth, while in the rear flank (left-hand side of the storm in Fig. 5a, which is the west) a turbulent cloud base could be seen, and the cloud base was rough. In the center of the storm, the cloud base extended lower towards the ground, and a funnel cloud with a diameter of 100 m touched the ground (Figs. 5a, d). It was accompanied by threedimensional streamlines swirling around a vortex line located in the center of the funnel cloud (Fig. 5d), indicating cyclonically rotating updrafts around the main body of the tornado.
The tornado structure was well-organized at t = 105 min, during its maturity. The three-dimensional pressure perturbation of the storm showed a clear column-like pattern. Outside the main part of the low-pressure center, in the mid-levels, a column of significant pressure drop was formed and connected to the ground (Fig. 5b). The inner  core of this pressure column showed a helical pattern, with stronger values at lower levels (Figs. 5b, e). The strongest pressure perturbation inside the funnel cloud exceeded 100 hPa. Another feature was the columns of vertical vorticity (Fig. 5c). Clusters of weaker columns existed at the rear flank, consistent with the occurrence of the turbulent cloud base (Fig. 5a). The location of significant vertical vorticity columns coincided with the column of strong pressure perturbation. Approximately two main columns of vorticity could be distinguished-one swirling around the other, and matching well with the helical pattern of the pressure field. The two columns of vorticity formed the double-helix structure of the vortex, which is further discussed in the following section. Similar to the pressure perturbation column, the vertical vorticity also had maxima close to the ground, indicating stronger rotation in the lower part of the funnel cloud (Fig. 5f).
The funnel cloud was surrounded by strong updrafts, with a maximum intensity at the boundary of the funnel  (d) is the threedimensional streamline (vortex line) that passes the rotation center at 100 m AGL. Panels (g)-(i) are two-dimensional contours of vertical velocity (w), horizontal velocity (U), and meridional velocity (v), respectively, on the plane denoted by the red dashed rectangles in panels (a)-(c). The wind vectors on the same plane are also shown by black vectors in panel (i). Please note that the isosurface in panels (a) and (d) are opaque, and the variation in colors only indicates the shades from camera light of the same value; while the isosurfaces in panels (b), (c), (d), and (e) are transparent; hence, the overlapping of isosurfaces may result in variation in colors but does not indicate a different value as well.
APRIL 2018 cloud, and decreasing both inwardly and outwardly (Fig. 5c). It also went down close to the ground. Two downdraft maxima could be seen inside the column-one between 600 and 700 m AGL, and the other close to the ground below 200 m AGL, with a diameter within the grid spacing of 25 m.
The full horizontal wind speed showed a similar pattern of maxima surrounding the funnel cloud, but with a maximum intensity close to the ground (Fig. 5h), as hypothesized according to Rankine vortex theory (Davies-Jones, 2015). This strong wind at the ground exceeded the EF4 threshold (74 m s -1 ). The meridional wind plot showed cyclonic rotation within the funnel cloud (Fig. 5i).
The wind field of the simulated tornado below the cloud base bore a strong resemblance to a tiny tropical cyclone, both in its rotational strong wind around the funnel cloud and the compensating downdraft in the middle. However, contrary to what one may have imagined, the funnel cloud was not accompanied by reflectivity echoes (data not shown), since the funnel cloud consisted of cloud droplets (and the debris, which was excluded in this simulation), instead of hydrometeors with larger sizes, such as rain droplets, or others that may have reflected significant echoes to the radar.
Judging from the structures at various scales, it should be pointed out that the vorticity and wind field patterns associated with the visible funnel cloud of the tornado under the cloud base were merely a small part of the entire structure of the tornadic circulation. Adding even more complexity to the situation, the structure of the tornado was highly spiral-like (helical). This shows the complexity of tornado structure, and the need to illustrate the three-dimensional structure through isosurface analyses.

Double-helix structure of the vortex
An interesting finding regarding the simulated tornado was its double-helix structure in terms of the full, three-dimensional vorticity field. At t = 91 min, far before tornadogenesis, a quasi-horizontal branch (QHB) at the lower part of the vorticity column could be seen (Fig. 6a). Several minutes later, another quasi-vertical branch (QVB) of the vorticity column took shape, which then formed a double-helix structure, together with the QHB (Fig. 6b). After the formation of the double-helix Fig. 6. The double-helix structure of the tornado indicated by full, three-dimensional vorticity and streamlines at t = 91, 99, 105, and 111 min. structure, the QVB grew stronger (Figs. 6c, d), and was associated with the significant vertical vorticity column (Fig. 5c) and the funnel cloud (Fig. 5a) of the simulated tornado. The QHB gradually weakened but persisted (Fig. 6d) until the demise of the tornado, associated with the downward growth of the cloud base before the formation of the funnel cloud.
Examination of the three-dimensional streamlines may help in understanding the intensification of the QVB of the double-helix structure. At t = 99 min, a streamline that passed through the QVB that rotated around its own axis near the surface horizontally could be seen (Fig. 6b). This rotation was forced by the horizontal temperature gradient formed by the cold pool in the rear flank downdraft (RFD) region, forming horizontal vorticity through baroclinic generation (Wurman et al., 2012). This source of horizontal vorticity was then tilted to be vertical, and stretched to be significantly stronger when entering the strong updraft of the tornado. At t = 105 min, more streamlines that passed through the QVB could be seen, with stronger rotation around its axis (Fig. 6c). At t = 111 min, when the QVB became much wider, the rotational streamlines that passed through it became weaker at the surface, indicating a reduction in the supplement of the vorticity into the tornado. For the streamlines at upper levels that swirled around the tornado, however, the rotation was still strong, with enlarged width. It should be noted that the streamline cannot be equally interpreted as the trajectory, but it does reveal the distribution of the wind field at a specific time.
Although several studies on fluid dynamics and laboratory simulations of tornadoes have discussed this structure (Rotunno, 2013), there has been little discussion on the spiral double-helix structure of vorticity, in either observations or numerical simulation studies. What was labeled as streamwise vorticity current (SVC) in Orf et al. (2017) may well have been showing similar structures to the QHB in our simulation. Their SVC, however, is generated along the baroclinic gradient within the forward flank, not the rear flank. Further research should be carried out on the relationship between the QHB and SVC, the possible explanation of the double-helix structure, and the impact or implications that it may have on the formation of tornadoes.

Formation of the tornado and hazardous winds 4.1 Formation of the simulated tornado
Having revealed the detailed structure of the simulated tornado in its mature stage, a further and arguably more intriguing question arises as follows: how did its formation happen? The three-dimensional isosurfaces of variables including pressure perturbation, vertical velocity (the updraft intensity), vertical vorticity (the rotation intensity), and cloud water mixing ratio (the funnel cloud) help to reveal this process.
During the formation stage of the supercell (stage 3 in Fig. 3), the updraft motion was generally strong in the mid-level of the storm at t = 80 min (Fig. 7a). A lowpressure center with perturbations of -5 and -10 hPa could be seen below 500 m AGL and close to the ground (Fig. 7e). It was accompanied by 20-m s -1 updrafts that extended below 2 km AGL from the main updraft of the mesocyclone at higher levels (Fig. 7a). A cluster of positive and negative vertical vorticity columns could be seen at the RFD region at this time-some of them close to the ground, and others around 1 km AGL (Fig. 7i). The positive vorticity columns were generally stronger, but the quantity of negative vorticity columns was larger. The cloud base of the supercell in the area where updrafts developed was about 800 m AGL (Fig. 7m).
Immediately before the sharp intensification of the mid-level rotation (stage 4 in Fig. 3a), with an intensified low-level updraft exceeding 40 m s -1 at t = 90 min (Fig. 7b), the pressure perturbation core grew larger and stronger, and the -10-hPa core touched the ground (Fig. 7f). In terms of vertical vorticity, the total number of columns dropped, while several long and thin cores remained (Fig. 7j). The cloud base in the region of the pressure perturbation center grew downward to about 500 m AGL (Fig. 7n).
Before tornadogenesis (stage 4 in Fig. 3b), the strong updraft exceeding 40 m s -1 formed a long column from the ground to about 6 km AGL at t = 100 min, connecting to the updraft of the mid-level mesocyclone (Fig. 7c). Accompanying it were columns of pressure perturbation and significant vertical vorticity columns. An intensified pressure perturbation core of over -25 hPa took shape, but did not touch the ground (Fig. 7g). Aside from the main part of the vertical vorticity column that accompanied the significant updraft, a thinner but stronger branch of vertical vorticity formed to its forward flank, which touched the ground, and was associated with the tornado that was about to form (Fig. 7k). The cloud base of the low-level mesocyclone grew lower to about 400 m, but the tornado funnel cloud did not form until later (Fig. 7o).
After tornadogenesis (stage 4 in Fig. 3b), the strong updraft extended close to the ground at t = 110 min (Fig. 7d), and the columns of pressure perturbation and vertical vorticity became both taller and stronger. The in- Fig. 7. Evolutions of (a-d) three-dimensional vertical velocity, (e-h) pressure perturbation, (i-l) vertical vorticity, and (m-p) cloud water mixing ratio at t = 80, 90, 100, and 110 min. Please note that the overlapping of the transparent isosurfaces may result in variation in colors but does not indicate a different value. tensified -25-hPa pressure perturbation column touched the ground and extended upward to almost 4 km AGL (Fig. 7h). The vertical vorticity columns below the cloud base merged (Fig. 7l), and a funnel cloud with a diameter of 100 m formed (Fig. 7p). An interesting phenomenon was the formation of a series of vertical negative vorticity columns, which looked similar to the structure in the recent simulation study (Orf et al., 2017). Based on the location and pattern, these vorticity columns might have been formed by the shear instability at the boundary of the RFD cold pool and the warm environment. However, Orf et al. (2017) indicated that positive vorticity columns were also present and even outnumbered the negative vorticity columns (their Fig. 12). This difference awaits further investigation.
The trend of the related variables before and after tornadogenesis can also be seen in Fig. 8. The first strong signal in the vertical velocity field took place at t = 88 min, below 1 km AGL, near the cloud base (Fig. 8a), which then developed upwards. Two significant signals in the vertical velocity field could be found at t = 105-114 min and t = 128-136 min, relating to the first and second tornadogeneses, respectively. Most of the signals generally developed bottom-up with time, except for the strong near-surface updrafts at tornadogenesis, which were top-down (e.g., t = 103-104 min, Fig. 8a). The variation in pressure perturbation showed similar patterns (Fig. 8b), especially the clear signals near the cloud base and the general upward advection. One major difference between the vertical velocity and the pressure perturbation field was that the latter extended to the ground.
The vertical vorticity signals (Fig. 8c) appeared at the same time as the peak values of pressure perturbation (Fig. 8b). However, there were no early signals in the vertical vorticity starting at heights near the cloud base, as in the pressure perturbation field. In contrast, the vertical vorticity maximum emerged right next to the ground and developed upwards. This phenomenon agrees with previous observations of three tornado events in the United States, observed by using X-band Phased Array Mobile Weather Radar, in that "there is no evidence in any of the datasets that a TVS built downward from midlevels in the storm prior to or at the time of tornadogenesis" (French et al., 2013). Since radar observations can only sample specific elevation angles, our simulation, with dense layers below the cloud base, further supports their conclusion. After normalizing the vorticity values at each level (to reduce the difference in vorticity intensity between various heights), the signals of vorticity above the ground were revealed; however, they exhibited an upward advection pattern with time, with no obvious relationship to the near-surface rotations. In general, the trend of vorticity growth was bottom-up, with no clear sign of top-down growth.
Interestingly, the vertical velocity tended to correspond to the specific locations of maximum vertical vorticity (Fig. 8e). This field can be used to physically evalu-ate the intensity of the updraft that may stretch the strong vorticity, which is the general idea behind the calculation of updraft helicity at a specific level. The overlapping of strong updraft and rotation has also been noted in previous studies (Weisman and Rotunno, 2000). Both the upward and downward development of intense signals starting from the cloud base could be seen between t = 100 and t = 105 min, before tornadogenesis. This shows that, aside from the upward growth of updraft signals that intensified the mid-level rotation, tornadogenesis was triggered by the downward growth of strong vertical velocity impinging on the intense vorticity center. Meanwhile, after tornadogenesis at t = 105 min, although strong updraft centers were maintained close to the ground, the intense vorticity center was covered by strong downdrafts, which might be the reason why the tornado's structure was gradually broken. Similar trends were also found during the evolution of the second tornado (t = 128-136 min in Fig. 8e).

Formation of the hazardous winds
Besides the detailed structures of the tornado itself, another important insight of this work is how the hazardous winds were formed. There has long been a belief that wind damage during a tornado event is caused by the tornado, and a general principle of a tornado's onsite damage survey is to infer the wind intensity and tornado width from the severe damage and wreckage on the ground. However, our simulation showed that the cause of the hazardous winds might be more complex.
The simulated damage swath was approximately 30 km long and about 4 km wide (Fig. 9), lying from westsouthwest to east-northeast, which is quite similar to the results from the damage survey . The funnel cloud touched the ground and persisted for more than 10 min, leaving a much more severe damage path of about 9 km long and less than 1 km wide. The simulation results vividly showed the entire process of the formation and evolution of a funnel cloud (Fig. 9b). At t = 101 min 50 s, a cone-shaped cloud, in terms of a 1-g kg -1 isosurface of cloud vapor mixing ratio, started to form. It then moved downward rotationally and touched the ground at t = 104 min 20 s. During the following few minutes, the funnel cloud grew significantly wider, to a diameter of approximately 200 m. It became even wider at t = 110 min. The funnel cloud started to shrink at t = 112 min 30 s, and eventually vanished at t = 115 min (refer to the animation in the supplementary material).
It has long been known that the actual distribution of the damage exerted by a tornado will be asymmetric, even if the tornado circulation itself is symmetric about its center (Atkins et al., 2014). This is due to the forward motion of the tornado adding (subtracting) to the effective wind speed on the right (left) side, resulting in different ground-relative speeds. A corresponding conclusion from this hypothesis is that, in cases when the damage area is too large to be fully covered, or the right half is covered with too few damage indicators to be rated, as in this case, if the damage distribution of the tornado track on the left-hand side can be obtained, then the other half can be extrapolated without being overestimated, since that part should be equal or more severe. However, during the damage survey of the present case, we noticed that in some parts of the tornado track the left half of the motion was more severe than the right half, and the 4km-wide diameter of the measured tornado track was much larger than the tornado's circulation.
The simulation results provided a possible answer to this problem. As shown in Fig. 9c, there were three areas that exhibited significantly strong surface winds. Area 1 was the inflow area to the forward flank of the tornado. This area was accompanied by wind reaching the threshold of EF1, and was symmetric about the center of the tornado. Area 2 was the inflow area at the rear flank towards the tornado. This was collocated with the QHB of the double-helix structure of vorticity. The wind exerted by this area exceeded the EF2 threshold, and was located only in the left half of the tornado's track, which was asymmetric about the tornado center and was quite different from area 1. The third area was the main part of the tornado's circulation. The hypothesis of symmetricity caused by adding and subtracting of the rotational wind only applies to this region.
The authors are fully aware that the simulated surface wind speed may have been made artificially high by the use of the free-slip condition. Another uncertainty in this explanation comes from the difference in the low-level mesocyclone between the simulation and observations (Figs. 2b, e). This hypothesis regarding the three components of the surface hazardous wind requires further investigation.

Concluding remarks
Successful tornado simulations have long been a challenging problem in the realm of meteorological research, and none has been carried out for tornadoes that happened in China. This work has addressed this knowledge gap by using a numerical simulation performed in CM1. A strong supercell resembling the radar observations of the EF4 tornado that took place in Funing, Yancheng, China, was successfully simulated by using a WRF-based proximity sounding as the homogeneous background and the initiation of convection using the updraft-nudging method. The simulated storm produced two tornadoes with funnel clouds touching the ground and surface winds exceeding the EF4 threshold (74 m s -1 ), which coincided with the tornadoes in Funing and Sheyang County associated with the same event.
The detailed structures of fields including vertical vorticity, pressure perturbation, vertical motion, and streamlines accompanying a tornado funnel cloud, and their evolution during tornadogenesis, have been reported for the first time in the literature on tornadoes that occurred in China. Maximum surface vertical vorticity has been shown to be a good indicator of the tornado's evolution in the simulation, accompanied by peak values of negatively correlated trends in maximum surface wind velocity and minimum surface pressure. The widely used updraft helicity was shown to be representative of a midlevel mesocyclone instead of the tornado itself. The signal of tornadogenesis was found first at the cloud base in the pressure perturbation field, and developed both upward and downward in terms of maximum vertical velocity overlapped with the intense vertical vorticity centers. The maximum vertical vorticity was found to have emerged close to the ground and developed from the bottom up, with no obvious sign of top-down growth, consistent with previous observations by X-band Phased Array Mobile Weather Radar in the United States. The tornado's demise was found to accompany the strong downdraft overlapped with the intense vorticity centers. Also found was a double-helix structure of the full vorticity, which, to the best of our knowledge, had only been previously discussed in theoretical derivations of tornado fluid dynamics, and not captured in a real tornado case. A hypothesis is put forward that the hazardous winds in this tornado event may have been produced under the interaction of three distinct areas-more complex than the tornadic circulation itself.
One interesting finding of the present work is that a violent surface vortex, accompanied by intense radial inflow and a corner region, was able to be generated and maintained even though the simulation employed a freeslip lower boundary condition (the shear stress at the surface was set to the stress at the first interior grid point, as opposed to a zero-stress boundary condition). The success of this simulation (the simulated and observed storms had many striking similarities), despite the limited data and idealized simulation approach, is encouraging in that additional studies of historical tornadoes in China may be possible.
To overcome the limited representativeness of the initial environment, ongoing research aims for reasonable selection and comparison of the proximity soundings. Future work will include the influence of the initial proximity sounding on the development of tornadic-scale features, detailed analyses of tornadic-scale features that lead to tornadogenesis and failure, an explanation and application of the double-helix structure, and implementation of this method in more Chinese tornado cases.