Previous numerical simulations have shown that tornadoes are sensitive to the characteristics of near-surface inflow (Lewellen et al. 1997;Lewellen et al. 2000;Lewellen and Lewellen 2007a,b), which in turn is influenced by characteristics of surface terrain. While understanding of tornado behavior has improved over the last few decades, the impact of complex topography on tornado dynamics remains relatively unknown. Though field projects such as Verification of the Origins of Rotation in Tornadoes Experiment-Southeast (VORTEX-SE) have targeted tornado observations in regions of complex terrain, namely, the southeast United States, obtaining a comprehensive radar dataset of a tornado traversing a significant terrain feature presents several major challenges. In addition to beam tilt and Earth curvature preventing data collection near the surface, larger terrain features would likely produce beam blockage problems and increase the difficulty in determining radar sites. High-temporalresolution mobile radar data (Wurman et al. 1997;Pazmany et al. 2013;Kurdzo et al. 2017) would also be a necessity as terrain features (slope, height, etc.) may vary rapidly with respect to tornado location. Thus, the number of observational studies relating the effect of terrain features on tornado or supercell dynamics is limited, primarily restricted to inferred wind characteristics from analysis of environmental factors (e.g., Bosart et al. 2006;Schneider 2009;Katona et al. 2016) or damage surveys (e.g., Fujita 1989;Karstens et al. 2013).
Several case studies have examined how terrain may lead to more favorable conditions for the intensification of the supercell and mesocyclone. Schneider (2009) analyzed three separate tornado events that occurred over the Great Tennessee Valley and noted multiple effects from surface terrain that led to a higher likelihood of tornadogenesis, including increased low-level convergence, upslope flow resulting in a strengthened updraft and mesocyclone, pooling of low-level moisture yielding low LCL heights, backing of surface wind through valley channeling resulting in larger storm-relative helicity values, and forced stretching of vertical vorticity when moving from higher to lower terrain. Bosart et al. (2006) also analyzed a tornadic supercell over Great Barrington, Massachusetts, and found that fluctuations in mesocyclone strength aligned with patterns in complex terrain, and hypothesized that flow channeling from topographic configurations increased tornadogenesis likelihood. Similarly, Tang et al. (2016) observed that flow channeling within the Mohawk Valley in New York led to a maximum in moisture flux convergence and a subsequent rapid strengthening of the supercell as it moved into the valley. Climatological studies have identified that rugged terrain may reduce tornado frequency (Gaffin and Parker 2006) while plateaus may increase low-level winds and tornadogenesis frequency (e.g., northeast Alabama; Lyza and Knupp 2018). Markowski and Dotzek (2011) analyzed thermodynamic environments of simulated supercells with and without terrain-airflow forced over or around terrain was found to induce thermodynamic heterogeneities in storm environments.
On the tornado scale, damage surveys have been utilized to analyze terrain effects on tornadic wind distributions. Fujita (1989) analyzed the damage swath of the 21 July 1987, Teton-Yellowstone, Wyoming, EF4 tornado and found the tornado weakened while ascending the 3000-m-high plateau. More recently, Karstens et al. (2013) analyzed tree fall patterns associated with both the 22 May 2011 Joplin, Missouri, EF5 tornado and the 27 April 2011 EF4 Tuscaloosa-Birmingham, Alabama, tornado. In the Tuscaloosa area, the study noted enhanced tree fall swaths to the left of the tornado far from the tornado's center and parallel to valley channels, which suggests that flow channeling of inflow led to increased tree fall. A shortcoming in studying terrain effects through damage surveys, however, is that variations in both the density of damage indicators and structural resistance makes it difficult to compare estimated winds along the track. Lewellen (2012, hereafter L12) presented a study of terrain effects on an idealized vortex using a large-eddy simulation (LES). The vortex was subjected to various terrain features including two adjacent hills, a valley oriented perpendicular to the translation of the vortex, a double ridge, a gapped ridge, and a slanted ridge. L12 noted surface terrain induced vortex tilt, path deviations, and overall variations in structure and strength of the vortex. For example, as the vortex ascended (descended) the ridge, the lower end of the vortex accelerated (decelerated) and tilted forward (backward) such that the central axis of the near-surface vortex was perpendicular to the terrain surface. Additionally, L12 noted a deflection of the vortex to the left (right) relative to vortex motion as the vortex ascended (descended) the ridge.
The present study aims to expand upon results presented in L12. To our knowledge, a comprehensive analysis of vortex sensitivity to both small and large variations in surface terrain has yet to be completed. The suite of terrain features in this study is broad, yet carefully chosen in order analyze variations in vortex behavior resulting from different types of terrain (e.g., hill vs valley), as well as subtle changes to similar terrain (e.g., small vs large hill). The overall goal is to address 1) what, if any, effects over different types of terrain can be generalized, 2) what terrain altitude changes are necessary to affect tornadoes, and 3) the overall variability in near-surface wind speeds and damage potential over varied terrain. It is important to point out that while the Great Plains is generally characterized by flat topography, subtle terrain features comparable to the ones presented in this study are common, while larger simulated terrain features might be found in other tornado-prone areas (e.g., southeast United States). The vortex and surface terrain are simulated through an LES with high spatiotemporal resolution. Lateral and upper boundary conditions are responsible for simulating a tornado-like vortex while terrain features are implemented through an immersed boundary method (IBM; Saiki and Biringen 1996), which relaxes all components of groundrelative flow at the terrain points to the local translation velocity.
The LES model is based on the computational fluid dynamic code created at the Research Institute for Applied Mechanics Computational Prediction of Airflow over Complex Terrain (RIAM-COMPACT; Uchida and Ohya 2003;Maruyama 2011) at Kyushu University. Specifications for the stretched meshed grid, time step, and model output are given in Table 1. Note that in the vertical direction, the grid spacing is constant at 2.5 m up to z 5 250 m to ensure sufficient capturing of terrain features before being stretched to 75.2 m at the top of the domain. Simulations with a coarser horizontal grid resolution of ;5 m have been performed with no appreciable change in qualitative results.
Flow within the open-domain LES is split into a grid scale (GS) and subgrid scale (SGS) through a filtering technique. The filtered GS mass continuity and Navier-Stokes equations (Maruyama 2011) are given by
where u i is the filtered wind component in the ith direction and F i represents the feedback force from the IBM, which relaxes ground-relative flow at the boundary to the local boundary velocity. Additionally, r is air density, p is filtered pressure, and S ij is the filtered deformation term. The eddy viscosity term n e is the summation of the viscosity and SGS viscosity. The LES boundary conditions, adapted from Bodine et al. (2016), generates a background vortex with the lower boundary translating in the negative x direction at a constant specified speed to simulate forward movement of the vortex. The dynamical lower boundary condition differs between the spinup period and the subsequent terrain period. First, the spinup period is run for 915 s and allows the vortex to reach a quasi-steady state before introducing terrain. During the spinup period, the semislip lower boundary condition imposes a surface roughness of z 0 5 0.001 m. After the spinup period, terrain is introduced from the right side of the domain and the lower boundary condition is defined through the IBM (discussed in further detail in the next section) rather than an imposed z 0 .
The lateral boundaries provide the source of angular momentum and inflow through a 200-m inflow layer of height h inf . Below h inf , nearly axisymmetric winds are imposed on the lateral boundary conditions to represent the swirling inflow at the mesocyclone scale. The wind profile is logarithmic up to h inf with a maximum value of 24.3 m s 21 in order to provide consistency with the imposed semislip boundary condition during the spinup period (Fig. 1a). The resulting mean angular momentum over the 200-m depth is ;10 000 m 2 s 21 . Above h inf , the lateral boundary conditions are imposed as no-slip with no flow entering or exiting the domain. The upper boundary condition induces an axisymmetric updraft with a radius of ;450 m and vertical velocity of 22.1 m s 21 that is responsible for driving flow from the lateral boundaries inward and upward (Fig. 1b). Placed within the central updraft is an imposed axisymmetric central downdraft; from R 5 172 m (one-third of the way from the center to edge of the domain) to 0 m, the vertical velocity decreases linearly from 22.1 to 220 m s 21 , respectively. While the imposed downdraft is not a necessary condition, it does aid in producing a vortex with more intense and numerous subvortices. Flow configurations are such that the mass flux leaving the domain through the updraft and the mass flux entering the domain through the inflow-transition layer are equal; i.e., the model imposes a net zero mass flux. In the cases where the inflow layer is partially blocked by terrain, the inflow layer is adjusted upward accordingly such that the mass flux into the domain still corresponds to that leaving the domain.
It should be noted that properties of the imposed wind profile on the lateral boundaries such as a deep log layer profile and an abrupt cutoff of flow at 200 m AGL, as well as the lack of adjustment to vortex motion and terrain intersection, may be different than what may physically develop from a strong vortex over a surface with lower roughness. Additionally, limited domain size and close proximity of the lateral boundaries to the vortex restricts the opportunity for the flow off lateral boundaries to adjust to the surface terrain before reaching the vortex. Nevertheless, vortex behavior within the parameter space in this study yields notable results that are worth documenting. Modifying wind profiles on the lateral and upper boundary conditions and increasing the domain size to allow for increased flow modification is left for future work.
During the terrain period, the IBM defines the lower boundary condition, details of which can be found in Goldstein et al. (1993) and Saiki and Biringen (1996), and are summarized below. A higher-resolution ''terrain grid'' (1000 3 1000 grid points) is defined such that sharp changes in the terrain feature can be adequately represented and transitions as the terrain translates across the model domain are smooth. The feedback force F i imposed on the terrain grid points x s can be represented mathematically as
where U i is the fluid velocity and y i is the specified constant translational velocity of the lower boundary. The coefficients a and b are dampening coefficients which govern the magnitude of numerical oscillations produced by the addition of the feedback function (Goldstein et al. 1993). In Saiki and Biringen (1996), a is chosen to be 24000 and b is chosen to be 260. After a range of larger and smaller values were tested, it was found that the best convergence and stability occurred for similar values; thus, these same values are used in the present study. The feedback force for each LES grid point is computed by linearly interpolating the feedback force at adjacent terrain grid points. Though the flow field in the near vicinity of the boundary contains inaccuracies as a result of first-order accuracy of interpolation, larger-scale flow is sufficiently captured for the purposes of the LES.
The terrain height at each grid point and time is given by
where x gp and y gp are the x and y position of each grid point, respectively, and f x and g y are the functions responsible for variation of terrain heights in the x and y directions, respectively. The terrain height for the 50-m ''base'' simulations for the four terrain categories-2D sinusoidal hill (hereafter sinusoid), 3D hill (hereafter hill), valley, and ridge-are shown in Fig. 2 with a base vortex translational velocity of 10 m s 21 . Unless otherwise noted, the terrain in the sinusoid and hill simulations repeat every one domain length. In total, a suite of 29 different simulations was run to evaluate different terrain-related factors that impact vortex structure (Table 2). Modifications from the 50-m base simulations include varying the maximum height of the terrain feature, changing the vortex translational velocity, steepening the approaching slope of the hill, offsetting the hill in the y direction, decreasing the frequency of successive hills, and oscillating the valley or ridge perpendicular to the direction of translation. Additionally, a control simulation, in which a zero-height terrain is defined, provides information on the background vortex produced by the boundary condition with no terrain influence. For brevity, horizontal and vertical winds are referred to as uy and w, with uy 10 and w 10 corresponding to the winds at 10 m AGL-components of the wind field are with respect to the z 5 0-m plane and quasi vortex relative (assumes the forward speed of the vortex matches the backward speed of the lower boundary). With regard to the former, it was found that while adjusting the vertical velocity to account for the gradual slope of the terrain does quantitatively affect the magnitudes of vertical velocity in some regions, especially for simulations where the terrain is steepest, qualitative results of w 10 , particularly with regards to areas of the most enhanced vertical motions, remained relatively unaffected. Since we felt that the vertical motions owing to terrain-following flow were important to incorporate in overall statistical calculations (section 4), we chose to stick with wind vectors that follow Cartesian coordinates. The track of uy 10 shows that the nearsurface background vortex has a relatively constant core radius of approximately 100 m, or a core diameter (D core ) of 200 m (Fig. 3). Relating uy 10 back to the EF scale, the background vortex corresponds to an upper-end EF-3, with clear swirl patterns of uy 10 . 60 m s 21 and small pockets of uy 10 . 70 m s 21 associated with secondary vortices (or subvortices) rotating about the main circulation, collocated with small areas w 10 . 25 m s 21 . The subvortices in the control simulation are quasi steady, with three to five separate subvortices existing at a given time, detectable through 3D visualizations of the pressure perturbation field (Fig. 4), calculated by subtracting the maximum pressure within the upper core region (R # 250 m and 750 # Z # 1250 m) from the pressure field at each analysis time. It is important to note that the perturbation pressure in Fig. 4 does not include the hydrostatic component of the true pressure field.
Last, the core diameter (D core ) of the background vortex can be a useful parameter to nondimensionalize spatial lengths of terrain features and domain size. When nondimensionalized by D core , the domain length in the x, y, and z directions is 5, 5, and 10, respectively. The terrain in the hill and sinusoid cases repeat every 5, and the diameter of the hill and width of the ridge (z terrain . 10 m) is just under 4. Last, the width of the valley (z terrain , 10 m) is approximately 1.5, just greater than the background vortex core diameter.
This section aims to highlight key vortex behavior from specific simulations within the hill, sinusoid, and valley sets. The chosen simulations for each set are as follows:
d Hill: hill_25m_base (hereafter H25), hill_50m_base (H50), hill_100m_base (H100), hill_50m_v20 (H50V), and hill_50m_steep2 (H50S)
d Sinusoid: sinusoid_25m_base (S25), sinusoid_50m_base (S50), and sinusoid_100m_base (S100)
d Valley: valley_50m_base (V50), valley_100m_base (V100), and valley_100m_snake100 (V100S)
Differences in uy 10 and w 10 distributions between the chosen simulations give insight on how vortex behavior changes with systematic modifications in terrain. A majority of the analysis is focused on 10 m AGL winds since estimated wind speeds corresponding to the EF scale are assessed at 10 m AGL. We have chosen not to highlight results from the ridge set in this analysis as the vortex displays the most notable behavior in the cases listed above. For discussion purposes, the terrain for the hill and sinusoid sets are split into three terrain segments: 1) the uphill segment, 2) the downhill segment, and 3) the gap segment located in between the downhill and uphill segments (z terrain , 10 m). Vortex quadrants relative to its forward motion-rear, front, left, and right-are also referenced in the analysis. Last, Lewellen et al. (2000) defined a corner flow swirl ratio S c , which is a dimensionless shape variable that is helpful in diagnosing low-level vortex behavior and is a function of the amount of low angular momentum flux into the vortex.
Decreasing S c from a high-swirl value intensifies the low-level vortex until its critical value, at which point further decrease results in a weakening low-level vortex.
The 10-m track analysis for H25, H50, and H100 highlights the influence of increasing terrain height on distributions of uy 10 and w 10 (Figs. 5 and 6, respectively). For H25, while values of w 10 are slightly higher in regions where z terrain . 10 m (Fig. 6a), there is an absence of notable patterns in uy 10 with respect to terrain (Fig. 5a). On the other hand, the vortex in H50 and H100 displays consistent regions of enhanced uy 10 located in the gap and uphill segments (see pink circles in Figs. 5b,c). Lower end uy 10 contours (30 or 40 m s 21 ) in Figs. 5b and 5c indicate that the vortex at 10 m AGL expands as it traverses uphill and contracts on the downhill, especially notable in H100. Additionally, areas of w 10 . 25 m s 21 are located almost exclusively on the uphill for H50 and H100 (Figs. 6b,c). Again, it is noted that while the decision to retain vertical velocity with respect to Cartesian coordinates rather than terrain-following coordinates (with w calculated as the velocity component normal to the local terrain surface) may be contributing somewhat to enhanced vertical motions here, areas of strongest w 10 still qualitatively mirror the results shown in Fig. 6; a comparison of w 10 using Cartesian versus terrainfollowing coordinates for H100 is shown in Fig. 7. Thus, while flow along terrain does contribute to some of the vertical TABLE 2. All simulations that are used for analysis sorted by terrain category. For brevity, terrain type is not included in the simulation name after the first row. The second term ''Xm'' after the first underscore (_) refers to the maximum height of the hill or ridge, or maximum depth of the valley. The last term after the second underscore refers to any changes made from the base simulation. Variations are as follows: ''vX'' 5 vortex translation speed is changed to X m s 21 , ''steepX'' 5 hill is steeper in the x direction with ''steep2'' being steeper than ''steep1,'' ''offsetX'' 5 center of hill is offset to the north at y 5 X m, ''spread'' 5 frequency of the hill is every two domain lengths as opposed to one, ''snakeX'' 5 valley or ridge oscillates in the direction of translation with a wavelength of one domain length and an amplitude of X m, ''curve'' 5 maximum depth of valley curves off out of the domain about halfway through the simulation.
Sinusoid set Valley set Ridge set velocity magnitude seen in Fig. 6, we remain confident that these areas of stronger w 10 are ''true'' enhancements, i.e., not just owing to upslope terrain.
Intense uy 10 and w 10 on the uphill segment are characterized by a noncontinuous swath, and are a result of transient, but robust, subvortices that are prominent on the uphill and suppressed on the downhill. Figures 8 and9 illustrate typical subvortex behavior during ascent. The preferred location of formation is in the rear quadrant in an area of relatively strong radial wind shears (Fig. 8a) above a broad area of 210-hPa pressure perturbation (Fig. 9a)-this is consistent with past studies which found that shearing instabilities play an essential role in subvortex formation (Rotunno 1978;Staley and Gall 1979;Walko and Gall 1984;Lewellen et al. 2000;Nolan 2012). As the subvortex rotates cyclonically around the main vortex, the strengthening subvortex (subvortex 2 in Figs. 8b,c) is characterized by 230-hPa pressure perturbation at 10 m AGL (Figs. 9b,c) with a component of motion radially inward. Additionally, both subvortices 1 and 2 at peak strength have a ''tail'' behind them that is associated with relatively large values of vertical vorticity generation both by stretching and tilting along areas of strong radial shears which could be vorticity sources for the subvortex. The subvortex terminates soon after being displaced from its initial radius (Fig. 8d), and the 230-hPa pressure perturbation isosurface lifts to 30-40 m AGL (Fig. 9d). Weakening of subvortices is rather abrupt, with uy 10 associated with the subvortex decreasing by as much as 30 m s 21 between subsequent 1.2-s model outputs.
Enhancement of w 10 on the uphill on the right quadrant of the vortex (Fig. 6) associated with terrain may enhance tilting and stretching of vorticity associated with subvortices-at the very least, shearing instabilities associated with radial gradients of vertical velocity are amplified in this region. The behavior of subvortices here is consistent with findings from Forbes (1998), who found a tendency for a tornado to produce suction vortices during an uphill climb, and Bluestein et al. (2018), who found that subvortices in the 31 May 2013 EF3 El Reno, Oklahoma, tornado typically formed in the left-rear quadrant and dissipated in the front quadrant, though that case occurred over relatively flat terrain. Compared to the control simulation, subvortices during vortex ascent are more robust, both in terms of strength (deeper pressure perturbations and stronger uy 10 /w 10 ) and size (larger volumes of negative pressure perturbations). We can also understand the physical mechanisms responsible for vortex contraction/expansion by investigating how distributions of ground-relative flow (i.e., low-and high-swirl regions) affect S c in Fig. 10 ( Lewellen et al. 2000;Lewellen and Lewellen 2007b;Lewellen 2012). As the vortex begins to ascend the hill, there is a pooling of low-swirl fluid on the lee side of the hill, which the vortex ingests as it crests the hill indicated by depleted ground-relative 10 m AGL wind vectors in Figs. 10a and10b. This ingestion of low-swirl fluid at the surface rapidly lowers S c and results in both an expansion of the lowlevel core radius and a one-cell vortex structure with stronger downdrafts restricted to farther off the surface (Fig. 11a). It is once again recognized that the proximity of the eastern lateral boundary to the vortex and its interaction with the introduced hill may be directly responsible for this pooling of low-swirl fluid; thus, it would be worthwhile to expand the domain and investigate modifications of areas of low-and high-swirl fluid.
Nevertheless, as the vortex descends the leeward side of the hill, the near-stagnant ground-relative flow ahead of the vortex is replaced by higher-swirl fluid, which rapidly increases S c of the vortex, leading to both contraction on the downhill (Figs. 10e,f) and transition to a two-cell structure where stronger downdrafts (w , 220 m s 21 ) impinges onto the surface (Fig. 11b). Additionally, there is a large area of w , 210 m s 21 that extends more than 300 m behind the vortex during the gap segment. This is in contrast to the control simulation, in which larger areas of w , 220 m s 21 generally remain above ;100 m AGL. The ingestion of high S c air results in repetitive areas of enhanced uy 10 when the vortex is in the gap segment.
There is also a tendency for the vortex in H50 and H100 to 1) deflect to the left (right) during ascent (descent) relative to its forward motion and 2) for the axis of the vortex to remain perpendicular to the surface, i.e., speed up and tilt forward (slow down and tilt backward) during ascent (descent), both of which are consistent with L12. Vortex tilt results in the central axis to be oriented perpendicular to the terrain surface, with more significant tilting generally restricted to below 600 m AGL (Fig. 12). Additionally, Figs. 12d-f also show descent of stronger (220-hPa) pressure perturbations from aloft down to the surface as the vortex descends the hill which is forced by the increase in S c and transition to a two-celled vortex during descent. Both deflections and tilting/accelerations are more significant for taller terrain (e.g., H100). L12 attributed the tilting of the vortex to the fact that the normal (to the boundary) pressure gradient force tends to zero in the Navier-Stokes equation, such that the central axis attempts to remain perpendicular to the terrain surface. Additionally, if the vortex is imagined as a line vortex, then impermeability at the surface dictates that there must be an image vortex below the surface such that relative angle between the vortex and the image vortex induces accelerations and deflections given by the Biot-Savart law.
Increasing vortex translation speed to 20 m s 21 (H50V) increases the magnitude of vortex contraction and expansion at 10 m AGL (Fig. 13b). Additionally, though the gap segment still shows a local maximum in uy 10 , the most notable area of enhanced uy 10 is positioned in first half of uphill segment (Fig. 13b) and is collocated with stronger w 10 (Fig. 14b). The faster translational velocity of the vortex may force higher swirl air ahead of the vortex on the uphill side to be deflected toward the vortex, increasing S c during the beginning of ascent before ingesting the lower-swirl air from the lee side of the hill. Additionally, because the central updraft is forced to move horizontally at 20 m s 21 , the lower-level vortex in H50V has more difficulty ''keeping up,'' and thus the magnitude of backward tilt is much more pronounced than the forward tilt. In fact, the vortex in H50V still maintains a slight backward tilt even during the uphill segment.
In H50S, which demonstrates the effects of a steeper approaching slope and a wider gap segment, largest uy 10 remains in the gap segment (Fig. 13c). Specifically, the largest area of enhanced uy 10 is located in an elongated eastwest band in the western half of the gap segment south of y 5 0 m, along with a smaller uy 10 maximum .70 m s 21 that occurs just before the uphill segment north of the vortex center at about y 5 100 m (Fig. 13c). Though the terrain is steeper, it is narrower in the x direction-thus, as the vortex enters the gap segment, the vortex has ample time to draw in high-swirl air from the lateral boundaries that is unaffected by terrain. This high-swirl air wraps around the vortex and results in the consistent swath of uy 10 . 60 m s 21 to the south of the vortex. Additionally, strong w remains restricted to the uphill portion of the terrain, perhaps even more so in FIG. 10. Sequence of 3D images from the hill_100m_base simulation showing vortex evolution over one hill from t 5 165.8 to 261.5 s. Shaded isosurfaces and arrows are as in Fig. 4, except that the depicted wind vectors are ground relative rather than vortex relative to highlight areas of low-and high-swirl fluid. Terrain height shading is as in Fig. 9.
H50S than in H50 due to greater steepness of the slope (Fig. 14c).
The simulations chosen in the sinusoid set for detailed 10-m analysis are S25, S50, and S100. Recall that for sinusoid simulations, terrain height does not vary in the x direction and extends to the lateral boundaries. Thus, comparisons of sinusoid to hill simulations aims to address the impact of increased areal extent and two-dimensionality of the terrain feature on vortex behavior. The uy 10 track has many of the same features seen in the hill simulations, including large swaths of enhanced uy 10 in the gap segment, as well as a broken pattern of stronger uy 10 coupled with higher w 10 from strong subvortices on the uphill segment (Figs. 15 and16). However, whereas the uy 10 track for H25 lacked any notable patterns, S25 does exhibit slightly larger swaths of uy 10 . 70 m s 21 in the gap segment, highlighting increased terrain influence on uy 10 patterns for sinusoid simulations. Last, as maximum terrain height increases, w 10 becomes more intense and confined to the uphill segment (Fig. 16).
The vortex in the sinusoid simulations is characterized by an increase in the variability of low-level vortex structure, including more significant expansion and contraction, larger deflections to the right and left of forward motion, and increased tilt and lag with respect to the central updraft. Figure 17 quantifies the lag of the vortex in hill_50m_base, hill_50m_v20, and sinusoid_100m_base by tracking the x position of the center of the vortex relative to the domain center. The center of the 10-m vortex for each time is defined as the location of the minimum in the averaged pressure perturbation field-the pressure perturbation field must be averaged spatially and temporally using the nearest neighbor average to capture the primary 10-m vortex position and eliminate influence of subvortices.
Figure 17 confirms that in H50 and S100, the vortex at 10 m decelerates (accelerates) while descending (ascending) the terrain, consistent with L12, with the lag much more pronounced in S100-maximum lag in H50 is about 40 m while lag for S100 is as much as 80 m. In H50V, deceleration (or lag) occurs from the top of the hill to the middle of the downhill segment and from the base of the hill to the middle of the uphill segment, while acceleration occurs elsewhere. As stated previously, the vortex lags the center of the domain nearly the entire H50V simulation, as the low-level vortex struggles to realign with a faster ''mesocyclone.'' In the majority of simulations, the 10-m vortex acceleration is sufficient such that it realigns with the central updraft. However, in simulations where the lag is most pronounced (e.g., S100 and H50V), the forward acceleration during descent is insufficient for realignment.
Figure 18 depicts the 10-m AGL vortex-relative winds as the vortex descends the terrain in S100 with the accompanying 3D visualization shown in Fig. 19. During descent, the vortex ingests high-swirl air ahead of it, increasing S c , rapidly contracting (Figs. 18a,b), and deepening near-surface pressure perturbations (Figs. 19a,b), similar to the process seen in the hill_100m_base simulation. Note that as a result of the terrain extension to the boundary, east-to-west (west-to-east) flow to the south (north) of the vortex has no option but to ascend the terrain, leaving behind pools of near-surface low-swirl fluid both to the northwest and southeast of the vortex, clearly discernible in the wind vectors in Fig. 19a. Ingestion of the low-swirl fluid may lead to the reexpansion of the vortex seen from Figs. 18b to 18c. As the circulation lags behind the center of the domain, it is unable to efficiently draw in the high-swirl air ahead of it which slowly depletes the angular momentum within the vortex (Figs. 18c and19c). Eventually, the old circulation is embedded in weaker flow and a new, more central circulation surrounds the old circulation (Fig. 18d). While the old circulations is actually associated with larger pressure perturbations in Fig. 19d, this quickly weakens within the next time step (not shown). It is important to note that the reorganization of the vortex within the center of the domain is almost certainly influenced by the lateral boundary conditions which supply a constant source of angular momentum toward the domain center. It is unclear whether or not the vortex would dissipate after being displaced so far from the ''mesocyclone'' without such boundary conditions.
The simulations within the valley set are V50, V100, and V100S (snaking valley). As the vortex enters the valley in V50 and V100, there is a slight deflection to the right (also noted in L12), with frequent north-south oriented enhanced uy 10 swaths within the valley associated with downslope flow off the southern side of the valley in the front-right quadrant (Figs. 20a,b). While uy 10 swaths greater than 70 m s 21 are restricted almost entirely to the south of y 5 0 m for V50, swaths in V100 extend slightly farther to the north. Additionally, there is a quasi-continuous stream of intense radial inflow from the rear quadrant of the vortex while in the valley. The tracks of w 10 for V50 and V100 reveal not only a dearth of longer-lived subvortices, but also a suppression of large vertical motions (Figs. 21a,b). When the FIG. 12. View of pressure perturbation isosurfaces from the south looking north depicting the evolution of vortex tilt and descent of pressure perturbations from aloft. Values of the pressure perturbation isosurfaces are as in Fig. 4.
vortex is in the valley, any component of radial inflow will act to enhance the downslope component of the flow, such that consistent areas of strong w 10 are suppressed; the lack of strong vertical motion remains even when using terrain-following coordinates (not shown).
For V100S, areas where uy 10 . 70 m s 21 occur in conjunction with z terrain . 10 m are limited, even when the oscillation of the valley brings the 10-m terrain height contour close to y 5 0 m (Fig. 20c). Thus, even when the domain and valley center are not collocated, areas of stronger uy 10 still favor lower elevations, i.e., the base of the valley. Additionally, areas of strong w 10 in V100S favor the region where the valley is at its southernmost point, with values of maximum w 10 (.50 m s 21 ) greater than those in V50 and V100 (Fig. 21c). When the valley oscillates, terrain influences lead to reduced radial inflow in the rear quadrant when at its northernmost point (Fig. 22a). On the other hand, when the vortex is at the southernmost point in the valley, the curve of the terrain is ideal for enhanced rear-quadrant radial inflow (Fig. 22b). The pressure response to the flow is much deeper perturbations at the southernmost point, which leads to larger vertical motions. vortex encounters the first hill in H50, H50V, and H50S, the vortex deviates to the south before it reaches the hill top, then shifts back northward as the vortex descends the first hill (Fig. 23a). For subsequent hills, the vortex remains displaced to the north from the base to the top of the hill, followed by a steep southward shift on the downhill just after the hill top, similar to numerical findings from L12 and Ahmed and Selvam (2015). This pattern repeats for the remainder of the simulation. Steepening the approaching slope (H50S) results in a slightly larger southern deviation on some of the downhill segments. Increasing vortex translational velocity to 20 m s 21 (H50V) results in larger variability of the vortex path, likely partially influenced by persistent backward tilt of the vortex even during ascent which forces reorganization of the 10-m vortex, similar to sinusoid_100m_base. For S25, S50, and S100, the vortex path is defined by a northern (southern) displacement as the vortex ascends (descends) the terrain (Fig. 23b). Beyond the first hill, phasing of the vortex displacements is overall similar to the hill simulations. The magnitude of the deviation is proportional to the maximum terrain height-the average southern displacement for S100 is ;80 m, a ;60% and a ;130% increase from the maximum southern deviations for S50 and S25, respectively. This is consistent with physical reasoning given earlier in regard to the Biot-Savart law, as a steeper terrain angle would yield a greater angle of the image vortex below the surface. Additionally, vortex displacement for the sinusoid simulations is slightly larger for a given maximum terrain height compared to the hill simulations. For example, maximum southern displacements for S50 and H50 are 50 and 35 m, respectively. This corresponds to an increased deviation of approximately 40% for S50, demonstrating that two-dimensionality of the terrain feature leads to larger vortex deflections. In both hill and sinusoid simulations, the magnitudes of displacements to the south exceed those to the north.
For valley simulations, the vortex path for V50, and more notably V100, have a slight preference to the southern half of the valley (Fig. 23c). Additionally, for V100S, the oscillation of the valley does indeed induce oscillations in the vortex path. However, the phase of the vortex and valley in V100S differed by about a quarter of a wavelength; the southernmost point in the vortex occurs about 250 m before the southernmost point in the valley.
This section aims to address generalized questions on the effects of terrain on vortex structure and dynamics through composite statistics (e.g., which terrain category has the largest uy 10 and w 10 ). The full suite of simulations in Table 2 is used in this analysis, and the data have been grouped by terrain category as well as terrain segment (uphill, downhill, gap) for hill and sinusoid cases. While vortex behavior with respect to terrain is not exactly periodic, and in some segments need not display periodic tendencies at all, overall periodic tendencies are such that averaging cases based on terrain location do yield some notable results.
First, the relationship between the height (radius) of maximum tangential velocity Z max (R max ) is examined as a function of terrain segment. The wind field at each height is decomposed into tangential and radial components; different center points at each height and time are calculated to ensure accuracy of axisymmetric calculations due to vortex tilt. Then, for each time, the tangential velocity is averaged into 5-m range and height bins; i.e., the grid points that correspond to a 5-m interval of range and height (e.g., between 50-55 m in range and 25-30 m in height AGL) are averaged. Next, the maximum tangential velocity for each height bin (across all range bins) is recorded, giving a profile of maximum tangential velocity for each time. This is repeated for all times, and the 95th percentile of all the profiles is calculated to get one representative profile for each simulation-this is repeated to obtain a tangential velocity profile with respect to R.
For hill simulations, the 95th-percentile tangential velocity height profile reveals maximum tangential velocity was slightly stronger (weaker) on the downhill (uphill; Fig. 24). Thus, while track analyses demonstrate weaker winds at 10 m AGL on the downhill slope, peak averaged tangential winds occur above 10 m AGL on the downhill. Additionally, while Z max does not differ by terrain segment (consistently in the 30-35-m height bin), analysis of R max in Fig. 24 confirms expansion and contraction patterns discussed previously. As the vortex crests the hill, it begins to ingest low-swirl fluid on the lee side of the hill which is eventually replaced by higher-swirl fluid as the vortex descends; thus, across all hill simulations, the vortex on the downhill has, on average, a smaller core radius than on the uphill, before reexpanding in the gap segment. The same general pattern occurs for sinusoid simulations (not shown), but with a few notable differences. First, for all terrain segments, Z max did vary between terrain segments, from 30-35 m in the gap segment to 40-45 m on the downhill; this provides evidence that two-dimensionality of the terrain feature resulted in larger variability of Z max perhaps due to increased FIG. 17. The zonal displacement of the center of the vortex relative to the center of the domain for hill_50m_base, hill_50m_v20, and sinusoid_100m_base. The dotted blue (red) lines represent the top (base) of the hill. frictional effects. Second, while R max gives the same pattern of vortex contraction/expansion, the core radius on the uphill and downhill terrain segments was 10-15 m wider in the sinusoid simulations compared to the hill simulations.
To calculate the 95th percentile of uy 10 and w 10 (hereafter uy p95 and w p95 ), data from the stretched grid are linearly interpolated to a regular 2 m 3 2 m grid such that each grid point receives equal representation for calculations. Then, uy p95 and w p95 are calculated from all grid points where R , 250 m over the entire simulated track at 10 m AGL, organized by both terrain category and segment (Table 3). All results discussed below undergo bootstrap significance testing to determine whether uy p95 and w p95 values between various simulations and terrain segments are significantly different. To determine significance between two sets, a random sample of 5000 points (from all grid points with R , 250 m through all times, totaling ;20 million points for one simulation not broken by terrain segments) is drawn for each of the two sets, and the 95th percentile of each random sample is then calculated. This is repeated 1000 times such that the result is two distributions of 1000 different 95th-percentile calculations from each of the two sets. The p value is simply the percentage of overlap between the two distributions. In the following analysis, if the p value is less than 0.2, 0.1, 0.05, and 0.01 (i.e., a confidence level of 80%, 90%, 95%, and 99%), then the 95th percentile between the two sets are said to be possibly, likely, significantly, and almost certainly different, respectively. If the p value is greater than 0.2, there is no statistical difference.
By terrain category (Table 3), overall uy p95 is almost certainly larger for the valley set compared all other simulations except the sinusoid simulation, in which there is no statistical difference. Overall uy p95 for the sinusoid simulations is also significantly larger than the control and ridge simulation and is almost certainly larger than the hill simulation. Last, uy p95 for the control and ridge simulation(s) are possibly larger than the hill simulations. Comparing hill and sinusoid sets directly, uy p95 for the sinusoid simulations is almost certainly larger (;3-4 m s 21 ) than the hill simulations through all three terrain segments, with uy p95 almost certainly the strongest (weakest) in the gap (downhill) segment, indicative of strongest winds occurring in lower terrain. Further evidence of stronger uy 10 favoring lower elevations is given in Fig. 25, with nearly 19% of all area within the lowest 10% of the terrain containing at least EF2 winds, far greater than the lowest 10%-20% of terrain (which contains 11.5% EF21 winds). Additionally, 2% of all area within the lowest 10% of the terrain also contains EF31 winds across all times which, from Fig. 25, is the largest of any 10% normalized terrain height bin, although this distinction is not as drastic as with the EF21 winds (green line in Fig. 25). It is also noteworthy to mention that there is a tendency for larger areas of EF21 winds and EF31 winds at the top of the terrain, influenced by more intense swirl spots or subvortices that rotate around the top as seen in hill_100m_base (see Fig. 10c). In fact, over 21% of the area of the highest 10% of terrain through all hill and sinusoid simulations have at least EF2 winds, although this may be skewed since the area of the top 10% of terrain is much lower compared to lower percentage bins, and will always occur in the center of the domain for hill simulations. Note that all area calculations are projections onto the z 5 0-m plane, and do not take into account small increases in area due to sloped terrain.
For vertical velocities, the ridge set has significantly larger w p95 than the hill set and almost certainly larger than all other simulations. Additionally, the hill set has a possibly larger w p95 than the sinusoid simulations and is almost certainly larger than the valley set and control simulation. Last, the sinusoid set is possibly larger than both the valley set and control simulation; w p95 for the valley set is not statistically different than the control simulation. Weaker w p95 for the valley set was already alluded to in the w 10 track analyses, resultant from radial inflow enhancing downslope flow and suppressing large w-the opposite is true for the ridge simulations where any radial inflow acts to enhance upward motion. The hill and sinusoid cases are a middle ground, where terrain-following wind vectors equate to both upslope and downslope flow depending on the location of the vortex relative to the terrain. However, the vortex in hill simulations brings upslope flow to a localized point when the vortex crests the hill, leading to a possibly larger w p95 . By terrain segment, the uphill segment for both hill and sinusoid sets almost certainly has
The following analysis quantifies changes in uy p95 and w p95 for terrain modifications within the hill and sinusoid sets (Table 4). The gap segment contains the largest uy p95 for all simulations shown except for hill_25m_base and hill_50m_v20. Recall for the hill_25m_base, uy 10 track analysis showed little difference between the three terrain segments, such that uy p95 differs less than 0.7 m s 21 due to the small amplitude of the terrain whereas for hill_50m_v20, the uphill contains the strongest uy p95 . For nearly all simulations, the downhill has the weakest uy p95 -in none of the simulations presented is the downhill the strongest segment.
For hill simulations, increasing the maximum height of the hill from 50 to 100 m results in a possible strengthening (weakening) of uy p95 in the gap (uphill/downhill) segment. Recall that the enhanced uy 10 in the gap segment was hypothesized to be a response to continued ingestion of higher swirl air ahead of the vortex during descent and in the gap segment. However, strengthening ceases when increasing the maximum hill height to 150 m-rather, this results in 1) an almost certain drop in uy p95 on the downhill segment and 2) a failure of uy p95 in the gap segment to increase further. A similar pattern emerges for the sinusoid simulations-when increasing terrain height from 25 to 50 m, there is an almost certain increase of uy p95 in the gap segment. Further increase to 100 m, however, results in a possible drop in uy p95 on the downhill and a nonsignificant change of uy p95 in the gap segment. Thus, while increasing terrain height can intensify the uy 10 in the gap segment, there is a critical height for which the terrain can become too disruptive, leading to extreme vortex tilting/displacement of the low-level vortex from the ''mesocyclone,'' the inability of the main circulation to continue to draw in higher angular momentum air, and forced reorganization near the center of the domain. This critical point occurs at a lower terrain height for the sinusoid set, since terrain effects are larger for sinusoid simulations compared to hill simulations.
To quantify risk posed by tornadoes over varying terrain, the percent area greater than EF1 (38.44 m s 21 ) and EF3 (60.79 m s 21 ) thresholds are computed (Table 5). For all hill and sinusoid simulations presented, the percentage of EF3 (EF1) winds are greater than (less than or equal to) the control simulation, which has a percent area of 0.7% (66.7%). Thus, larger areas of severe damage can be expected with the introduction of complex terrain while the overall area experiencing some degree of damage is reduced. However, within individual terrain segments, the risk can vary considerably.
For example, all simulations exhibited substantially higher risk of EF31 winds in the gap region compared to the control simulation. In contrast, some uphill and downhill segments have larger or smaller areas of EF31 winds. For EF11 winds, all simulations and associated segments experienced a reduction in area compared to the control. Thus, identifying the ''highest risk'' area depends on whether risk is defined as the region with the strongest wind speeds (most intense damage) or the largest areal extent of tornadic wind speeds (total damage area). If the former (latter) definition is used, the control simulation would be classified as a lowest (highest) risk simulation.
Composite upward vertical velocity statistics are quantified through the area of w 10 . 10 m s 21 (hereafter updraft area). Updraft area for the control simulation is 5.1%, which is relatively small when compared to nearly all simulations in Table 6 since terrain influences are nonexistent. Terrainfollowing flow within the terrain simulations are more subject to upslope components within the vortex, leading to the greater updraft area. Increasing terrain height, translational velocity, and steepness of the hill all lead to an increase in updraft area (Table 6) since these modifications act to increase the component of upslope flow during ascent. On the downhill, however, a faster translational velocity acts to decrease the upslope flow located in the left quadrant, leading to a reduction in updraft area in hill_50m_v20 compared to hill_50m_base.
This study uses an LES with a medium-swirl background vortex to investigate the effects of surface terrain, implemented through the IBM, on vortex structure and dynamics. A suite of 29 simulations is grouped into one of four categories-hill, sinusoid, valley, or ridge-in addition to a control simulation in which no terrain is introduced. Overall, the addition of surface terrain acts to increase the variability of along-track 10-m horizontal and vertical winds. This variability becomes more pronounced with an increase in terrain height, two-dimensionality of terrain (sinusoid versus hill), and increased translational velocity as these modifications act to enhance variability of the swirl magnitude of air surrounding the vortex. This is consistent with observations by Cannon et al. (2016) and Zenoble and Peterson (2017), who also noted an increase in variability in tree fall patterns over complex terrain compared to flat terrain. More specifically, the presence of terrain induces path deviations, vortex contraction/expansion, and vortex tilt not present in the control simulation, consistent with findings from L12. Expansion and contraction patterns are owing to changes in S c of the vortex. As the vortex crests the hill top, it begins to ingest low-swirl fluid that pools on the lee side of the mountain, which acts to decrease S c and expand the near-surface vortex.
As the vortex descends the lee side of the hill, low-swirl fluid ahead of the vortex is replaced by higher-swirl fluid, acting to increase S c , contract the near-surface vortex, and transition the vortex to a two-cell structure with stronger downdrafts and lower pressure perturbations impinging onto the surface; in contrast, stronger downdrafts in the control simulation remained far above the surface (above ;100 m AGL). In addition to expansion and contraction patterns, the central axis of the vortex remains perpendicular to the terrain surface, accelerating (decelerating) and deflecting left (right) of forward motion during ascent (descent). Last, the most intense winds tend to favor lower elevations, either in between hills or the lowest point in the valley, similar to findings from Cannon et al. (2016), with another maximum found to exist at the top of the hills (Selvam et al. 2015;Flynn and Islam 2019). For hill and sinusoid simulations, enhanced uy 10 in the gap segment is hypothesized to be a response to continued ingestion of higher-swirl fluid by the vortex in the gap segment. Subvortices are found to be stronger and larger as the vortex ascends the terrain especially compared to the control simulation, consistent with damage surveys by Forbes (1998). While quasi-steady subvortices in the control simulation are associated with 10-m horizontal and vertical wind speeds of 70 and 30 m s 21 , respectively, subvortices in the terrain simulations are significantly more robust with horizontal and vertical wind speeds up to 90 and 50 m s 21 , along with deeper, and volumetrically larger, negative pressure perturbations. Subvortices in terrain simulations are also shorter-lived as opposed to quasi-steady subvortices in the control simulation. Both the presence of subvortices and the forced ascent of intense terrain-following corner flow contributes to the strongest vertical velocities being located on the upslope portion of the terrain. Vertical velocities also increase with modifications that further enhance upward terrain-following flow, including increasing terrain height, steepness of the approaching slope, and translational velocity.
The simulations presented herein represent a small subset of possible terrain parameters with only one type of vortex flow being considered. Environmental heterogeneities within an FIG. 25. The percent area of EF21 (green) and EF31 (orange) winds for all hill and sinusoid simulations with R , 250 m binned by every 10% of normalized terrain height. For example, the second green dot is approximately 12% of the area with R , 250 m and terrain height between 10% and 20% of the maximum terrain height for all simulations and times that have winds greater than the EF2 threshold. observed near-tornado environment are much greater than simulations presented in this study. As stated at the end of section 2a, it is uncertain how much the unvarying, axisymmetric, and perhaps unnatural (for a surface with low roughness) lateral and upper boundary conditions, which provide all sources of vorticity and angular momentum, as well as the limited domain with respect to the vortex size, influence the evolution of the vortex as it interacts with the terrain. It could be expected, however, that increasing the terrain influence relative to the vortex-either by decreasing the size of the vortex or increasing the areal extent of the terrain-would result in an increased magnitude of vortex disruption. Last, during the spinup period, the lower boundary condition specifies a surface roughness of 0.001 m; once terrain is introduced, the IBM does not explicitly define a surface roughness parameter. It is speculated that the inclusion of a surface roughness parameter would likely yield more significant vortex responses to the terrain features (L12) and may change the character of the vortex altogether, and thus should be considered in future simulations.
Thermodynamically, influences of terrain on variables such as moisture and CAPE, and how these variables affect vortex behavior are not included within these simulations. Additionally, because the domain does not extend past 500 m ahead of the vortex, the simulations do not take into account influences of the terrain on flow patterns further downstream (relative to vortex translation), which could modify vorticity and angular momentum budgets and consequently, could change vortex behavior. Also, it is possible that the effects of pressure perturbations within the vortex core may impact the upper core, especially for higher terrain simulations. In such cases, it is caveated that the fixed upper boundary condition would not realistically capture the feedback. In future simulations, a nested circulation configuration, such as that used in Lewellen (2014), which terminates the vortex prior to the top of the domain, should be explored. In conclusion, while the suite of simulations presented herein provide insight into vortex behavior given these particular boundary conditions, we must be cautious in extending generalized results outside of this parameter space.
While a broad spectrum of surface terrain is introduced in this study, the breadth of topographic effects on tornadoes is quite large. The implementation of real terrain data from observed tornado tracks and simulating damage patterns (e.g., tree fall; Lombardo et al. 2015) would allow for a more direct comparison from LES to observed results, and may aid in identifying areas where improvements to the model terrain representation are needed. Moreover, using high-resolution supercell simulations as a method of nesting variable boundary conditions within the LES may be needed to account for storm evolution in real cases. Finally, intercomparisons among terrain models should be undertaken, including an implementation of LESs with terrainfollowing coordinates rather than a traditional Cartesian grid.
Downloaded from http://journals.ametsoc.org/jas/article-pdf/77/10/3277/5000626/jasd190321.pdf by guest on 19 September 2020
S A T R I O E T A L .