## Abstract

The structure and intensity of tornado-like vortices are examined using large-eddy simulations (LES) in an idealized framework. The analysis focuses on whether the simulated boundary layer contains resolved turbulent eddies, and whether most of the vertical component of turbulent momentum flux is resolved rather than parameterized. Initial conditions are first generated numerically using a “precursor simulation” with an axisymmetric model. A three-dimensional “baseline” LES is then integrated using these initial conditions plus random perturbations. With this baseline approach, the inner core of the simulated vortex clearly contains resolved turbulent eddies (as expected); however, the boundary layer inflow has very weak resolved turbulent eddies, and the subgrid model accounts for most of the vertical turbulent momentum flux (contrary to the design of these simulations). To overcome this problem, a second precursor simulation is conducted in which resolved turbulent fluctuations develop within a smaller, doubly periodic LES domain. Perturbation flow fields from this precursor LES are then “injected” into the large-domain LES at a specified radius. With this approach, the boundary layer inflow clearly contains resolved turbulent fluctuations, often organized as quasi-2D rolls, which persist into the inner core of the simulation; thus, the simulated tornado-like vortex *and* its inflowing boundary layer can be characterized as LES. When turbulence is injected, the inner-core vortex structure is always substantially different, the boundary layer inflow is typically deeper, and in most cases the maximum wind speeds are reduced compared to the baseline simulation.

## 1. Introduction

Numerical simulations can play a critical role in studies of the near-surface wind fields of tornadoes. In particular, numerical simulations can produce complete datasets in space and time, which allow for calculations of time-averaged wind speeds, short-term gusts, and their variations in space and time. Such information is especially important to the structural engineering community. In contrast, observational data within tornadoes are typically incomplete. The temporal evolution of flow features within tornadoes can be especially difficult to study with present-day observational techniques, leaving open questions about how wind speeds within tornadoes vary in time.

To study the detailed structure of the wind field in tornadoes, while maintaining reasonable computational cost, many numerical studies are based on idealized representations of “tornado-like vortices”—that is, vortices that have the same overall properties as tornadoes in nature, such as size and intensity, but are not connected to a parent thunderstorm. Several high-resolution (i.e., grid spacing less than 10 m) three-dimensional modeling studies have been published in recent years using this type of approach. For example, Lewellen et al. (1997, 2000), and Lewellen and Lewellen (2007a,b) studied the detailed structure of winds in tornado-like vortices in a series of numerical experiments. Their model domain used specified inflow conditions at the lateral boundaries and specified outflow conditions at the upper boundary. This setup is based partly on laboratory “vortex chamber” experiments, but also approximately represents the conditions surrounding the low-level mesocyclone in a supercell thunderstorm [see, e.g., Rotunno (2013), for a review]. Most 3D simulation studies use this approach, or some subtle variation thereof (e.g., Kuai et al. 2008; Maruyama 2011; Natarajan and Hangan 2012; Liu and Ishihara 2015).

Most of the studies cited above used the modeling technique known as large-eddy simulation (LES). In principle, LES uses sufficiently small grid resolution to represent the largest and most energetic features in a turbulent flow. LES requires a subgrid turbulence model, which accounts for the effects of unresolved turbulence on resolved-scale fields. Subgrid models for LES are often designed based on theoretical conditions [e.g., chapter 6 of Wyngaard (2010)]. However, the underlying assumption of LES is that resolved fluctuations contain most of the turbulent kinetic energy and turbulent momentum flux, and that subgrid turbulence models have only a small influence on results [except very near the surface, where turbulent structures becomes vanishingly small (e.g., Wyngaard 2010, p. 215) and, thus, are typically parameterized in atmospheric LES].

This issue—that is, whether the simulated flow actually contains turbulent eddies, and thus can truly be called “large-eddy” simulation—seems to receive little attention in published literature on tornado-like vortices. In all the LES studies cited above, the vortex core (defined roughly as the region near or inward of the radius of maximum wind speed) clearly contains turbulent fluctuations. However, the boundary layer inflow (i.e., the roughly 100-m-deep layer next to the surface that flows inward toward the vortex core) seems to *not* contain turbulent eddies in most of these studies, at least based on our interpretation of figures in those articles.

The underlying issue is not just related to semantics (i.e., what it means to call something a large-eddy simulation). The modeling components typically used in atmospheric LES, such as the aforementioned subgrid models, are typically formulated under the assumption that most of the turbulent flux is represented via the model’s governing equations (except very near the surface). Thus, these models probably will not work effectively at representing *all* the effects of turbulence in flows that do not contain any turbulent fluctuations. For a clear example, see Markowski and Bryan (2016).

For simulations of tornado-like vortices, the properties of the inflowing boundary layer are especially important. Many previous studies (e.g., Lewellen and Lewellen 2007a,b) (as well as simulations herein) demonstrate that the structure of the inner core of tornado-like vortices, and thus maximum wind speeds, can be related to properties of the boundary layer inflow. Our goal is to handle this boundary layer inflow, as well as the tornado’s core, using LES (except for a thin layer near the surface to be handled by a surface-layer parameterization). As demonstrated below, our methodology also allows for mature, fully developed, resolved turbulent fluctuations to travel into the tornado core, which we suspect has been missing from some previous numerical studies of tornado-like vortices.

In this article, we describe a unique approach to studying tornado-like vortices using LES. To ensure that most of the boundary layer contains resolved turbulent fluctuations, while maintaining reasonable computational cost, we find it necessary to insert turbulence at the edge of a high-resolution subdomain. To this end, a fully developed turbulence field is developed in a separate simulation and is then “injected” within the inflow to the high-resolution subdomain. The inclusion of resolved turbulence has a notable impact on the simulated boundary layer inflow, and alters the low-level structure of the tornado-like vortex in our simulations.

This article is the second in a series that examines wind fields in simulated tornado-like vortices. The first, Rotunno et al. (2016, hereafter RBND), explained the overall motivation of this work and highlighted the need to study the role of turbulence in the boundary layer inflow using 3D turbulence-resolving simulations. RBND also established the overall framework for the simulations in this article, which are based on the “Fiedler Chamber” approach in which an updraft forcing is specified in the middle of a rotating domain (Fiedler 1994; Fiedler 1995; Rotunno 2013). Although this approach differs from most previous LES studies of tornado-like vortices (cited above), the overall conclusions of this study concerning the role of turbulence on the boundary layer inflow likely applies to other modeling approaches.

## 2. Model description

The physical setup for this study is described and justified in RBND. Briefly, a dry and neutrally stratified fluid is simulated in a closed domain that rotates at a specified rate . A forcing term is included in the vertical momentum equation as a surrogate for a convective updraft in a supercell thunderstorm, following Fiedler (1995). Surface drag is applied on the lower boundary, which allows air to flow into the center of the domain near the surface (e.g., Rotunno 2013). At some region near the center of the domain, air converges and then rises upward from the surface, into the region of updraft forcing. A deep (7 km) damping layer is applied in the upper part of the domain to gradually eliminate features that move upward and away from the primary region of interest for this study, which is the lowest 1 km above the surface.

For all simulations we use Cloud Model 1 (CM1), a nonhydrostatic, time-dependent numerical model. An advantage of CM1 for this study is the ability to integrate different equation sets. Here we use a set of equations for which density is assumed to be constant and which permit acoustic waves, although the three-dimensional divergence is always negligibly small; hence, RBND refer to these equations as “effectively incompressible” (see their section 3a for more details). We chose this equation set for several reasons, such as the simple form of the equations for dynamical analysis, and the relative ease of incorporating the equations into the CM1 solver for distributed-memory supercomputers. Nevertheless, the neglect of compressibility effects herein may affect some details, such as the intensity of small-scale subvortices (i.e., “suction vortices,” Fujita 1970) where instantaneous wind speed can exceed 200 m s^{−1} in our simulations.

For a portion of this work we use an axisymmetric set of these equations (i.e., equations on a cylindrical coordinate system in which azimuthal variations are neglected). The goal with this setup is to develop an approximately steady solution which is then used to initialize the three-dimensional large-eddy simulations, thereby avoiding a costly “spinup” period.

The primary set of simulations uses a Cartesian grid with gradually varying grid spacing. The entire domain is illustrated in Fig. 1, which shows every 10th grid point as a blue dot. For all simulations herein (including the axisymmetric simulations) the vertical grid spacing is a constant 2.5 m for 1 km; above 1 km, increases gradually to a maximum of 248 m at the top of the domain ( km). The horizontal grid spacing is a constant 5 m in a 4 km 4 km region at the center of the domain (Fig. 1b); increases gradually to 220 m at the lateral boundaries of the domain. The overall domain has 1120 1120 512 grid points. The axisymmetric grid is identical to the right half of Fig. 1a, and has radial grid spacing of 5 m for radius km.

In our previous study (RBND) we used a constant coefficient of viscosity in the stress terms of the governing equations and a no-slip lower boundary condition. A primary goal of that study was to place previous studies into a broader context by explicitly controlling the Reynolds number. For the present study we are interested in details of turbulent flow structures in the boundary layer of tornado-like vortices, which requires different formulations for the stress terms and boundary conditions. To maintain reasonable computational resources we choose to conduct LES (see section 1). As explained above, subgrid turbulence models for LES of the atmosphere typically presuppose that the simulated flow contains turbulent fluctuations, that is, eddies (except very near the surface). However, because of the limited region of high resolution (Fig. 1), we cannot conduct LES over the entire domain. Therefore, for the relatively coarse-resolution portion of the domain we parameterize all turbulence effects using a simple scheme that has proven to work reasonably well for tropical cyclones (e.g., Bryan 2012) and shear-driven atmospheric boundary layers (e.g., Markowski and Bryan 2016). Details of the two subgrid models, and the transition between them, are provided in sections 3 and 4 below. The eddy injection scheme, described in section 5, turns out to be a crucial component of the transition between these two modeling frameworks in this study.

For the lower boundary condition in all simulations we use a traditional “semislip” type of boundary condition for atmospheric LES in which the surface stress is parameterized based on a logarithmic-layer equation. Specifically, surface stress is determined from the predicted horizontal wind speed at the lowest model level () and a specified roughness length :

where is the height of the lowest model level ( m herein), is the von Kármán constant, and is the friction velocity. It should be noted that turbulence near the surface in most atmospheric LES (including the simulations herein) is largely parameterized because turbulent eddies are smaller near rigid boundaries. For a review of the topic, and various methods that can be used to represent turbulence near surfaces, see Piomelli and Balaras (2002).

A free-slip condition is used for all lateral boundaries, and a no-slip condition is used for the upper boundary. Results are not sensitive to the boundary conditions at the top and/or sides of the domain, primarily because the damper layer (discussed above) is effective at removing flow features from upper levels, which are irrelevant to this study (see also the appendix in RBND).

Although the lower boundary is flat and uniform (i.e., is constant) for this article, these conditions are not required. By varying (in *x* and/or *y*), variable surface conditions could be investigated. CM1 also has terrain-following coordinates that could be used to study terrain effects. Different types of canopy models could, in principle, be implemented to more rigorously study forested, or perhaps urban, surface conditions. Nevertheless, we defer these topics to future work.

The primary solver of CM1 uses a three-step Runge–Kutta time integration scheme with a split-explicit method to account for acoustic modes (Wicker and Skamarock 2002). We use a fifth-order weighted essentially nonoscillatory (WENO) scheme (Shu 2001; Borges et al. 2008) for advection terms, and second-order centered differencing on a staggered C grid for all other derivatives. The time step is adaptive, meaning that changes throughout the simulation to maintain numerical stability. The scheme for this study requires the maximum Courant number to remain less than 1, and further checks that horizontal diffusion tendencies are numerically stable. The latter condition requires where *K* is an eddy viscosity and (which varies across the domain, e.g., Fig. 1). Vertical diffusion is calculated using an absolutely stable Crank–Nicholson scheme as described in Bryan (2012). In the following simulations, is as small as 0.01 s.

Finally, we note that airborne debris is not considered herein. Recent studies by Lewellen et al. (2008) and Bodine et al. (2016) present two different methods for incorporating debris effects into LES, which generally is expected to reduce near-surface wind speeds.

## 3. Precursor simulation 1: Axisymmetric model

Although it is possible to spin up a tornado-like vortex using the 3D model, we find this method to be prohibitively expensive. Instead, we integrate the axisymmetric version of CM1 until a strong and approximately steady vortex has developed. Because of the lower expense (the total cost of a 3D LES herein is reduced by roughly a factor of 7), this approach allows us to explore a broader range of the parameter space, including specified surface roughness length , swirl ratio , and forcing velocity scale *W* (see RBND for definitions of and *W*). Some results from sensitivity experiments, in which grid resolution and other model parameters are varied, are reported in a companion article (Nolan et al. 2017). All simulations herein use = 0.2 m, , and m s^{−1}.

### a. Model configuration

The governing equations of the axisymmetric model are listed in RBND, although here we replace the stress terms [next-to-last terms of Eqs. (1a)–(1c) of RBND] with the following parameterization for turbulence:

where and are the radial and tangential components of velocity, respectively; and and are the horizontal and vertical eddy viscosities, respectively. The two eddy viscosities are determined from separate relations:

where is a (constant) length scale for the radial direction, and is a length scale for the vertical direction determined from the following relation:

(e.g., Mason and Thompson 1992), where is a (constant) asymptotic length scale. The constants and have been shown to have a substantial effect on the structure and maximum intensity of axisymmetric tropical cyclones (e.g., Rotunno and Bryan 2012). The same sensitivity exists for axisymmetric tornado simulations, which we have confirmed by some limited sensitivity simulations (not shown). Herein, we use m, which produce reasonable structures in the low levels of our simulations, as shown below.

### b. Results

Time series of maximum vertical velocity and maximum tangential velocity from the axisymmetric simulation are shown in Fig. 2. For the settings used here, an approximately steady value of is reached early in the simulation (by 1000 s) whereas the maximum tangential velocity requires much longer (4000 s). Thereafter, the essentially steady character of the simulation is apparent in Fig. 2, although we note that the radius of maximum wind speed slowly increases throughout the entire simulation (e.g., from 197 to 217 m between 15 000 and 20 000 s) and, hence, the simulation is not exactly steady.

Figure 3 illustrates the simulated flow fields over the entire domain and within the inner-core region (inset). The height of is 95 m, which is approximately the same as the depth of the inflow layer. The radius of is 240 m. The inner-core vortex structure has a clear inertial oscillation with height attributable to inertia waves (e.g., Batchelor 1967, p. 559), which is a common feature of strongly rotating axisymmetric simulations (e.g., Rotunno 2014).

Time-averaged output from the axisymmetric model is used to initialize the 3D simulations that are the focus of our study. For the cases shown below, we average the three components of velocity and the pressure variable *ϕ* using output fields from the axisymmetric model written every 100 s for s.

## 4. Baseline simulation

### a. Model configuration

The three-dimensional (Cartesian coordinate) equations for this study (using tensor notation) are as follows:

where is the pressure variable (*p* is pressure and *ρ* is density), is rotation rate, *F* is the updraft forcing term (see RBND), is the stress tensor, is a relaxation time scale, is a parameter that gradually increases with height to smoothly turn on the damping layer (see RBND), and is the speed of sound. We use s and m s^{−1} herein.

The stress is determined in two ways, depending on location within the domain. For most of the inner fine-mesh (which is shown by a red box in Fig. 1) subgrid stress is parameterized in a standard manner for LES (e.g., Moeng and Sullivan 2015):

where *K* is an eddy viscosity determined from subgrid turbulence kinetic energy :

with . Following Deardorff (1980), is determined from a time-dependent equation:

where . The length scale *l* is limited near the surface:

where and .

Outside of the fine-mesh subdomain, the parameterization scheme is intended to account for all effects of turbulence (i.e., turbulent fluctuations are presumed *not* to exist), and is essentially the same as the scheme described in section 3. For completeness, the formulation is as follows:

where the eddy viscosities are determined from the relations:

The length scales and are identical to those in section 3.

Because an LES scheme is used for on part of the domain, and a scheme for completely parameterized turbulence is used for everywhere else, it was necessary to develop a method to smoothly transition between these two schemes. Figure 4 illustrates our method. Assuming the fine-mesh grid is 4 km 4 km, we specify a “transition zone” at a radius of 2 km from the domain center (i.e., a ring-shaped transition zone is placed inside the square-shaped fine-mesh subdomain). As illustrated in Fig. 4a, the LES subgrid model constant is reduced to zero on the *inner part* of this transition zone using a hyperbolic tangent function. The parameter is varied consistently. For the parameterized turbulence scheme, the length scales and are increased smoothly on the *outer part* of the transition zone (Fig. 4b) using a hyperbolic tangent function. Although all simulations in this article use m, we show an example with m in Fig. 4b to illustrate how these two variables need not be identical.

As a consequence of these formulations—and particularly the decision to vary one set of parameters on the inside of the transition zone, while another set of parameters vary on the outside of this zone—there is no subgrid turbulence model active in the middle of the transition zone. This setup is important because turbulent fluctuations (i.e., resolved eddies) must develop in the transition zone to replace the turbulence parameterization, (10) and (11), that is turned off in the same region. Because subgrid turbulence models are typically dissipative, it is advantageous to have essentially no turbulence model in the middle of the transition zone, thus minimizing dissipation and allowing perturbations to grow as rapidly as possible. A similar approach was used by Muñoz-Esparza et al. (2014, p. 418).

In contrast to most recent studies with embedded LES domains (e.g., Mirocha et al. 2014; Muñoz-Esparza et al. 2014; Munters et al. 2016), herein the outflow from the fine-mesh portion of the domain is in the *vertical* direction (rather than at a lateral boundary, or at a region of a domain where horizontal grid spacing varies). Consequently, variations in , , , and are made in the vertical direction in the same manner as the inflow “ring” illustrated in Fig. 4. For the model domain used herein (Fig. 1) the vertical transition zone is centered at km.

For initial conditions, the time-averaged results from the precursor 1 (axisymmetric) simulation are interpolated to the 3D Cartesian grid. Random perturbations are then added to initiate three-dimensional motions; specifically, random perturbations of 0.25 m s^{−1} are added to the horizontal velocity components at every grid point. (These perturbations do not satisfy incompressibility, and the *ϕ* field is not modified based on these perturbations; acoustic waves act to adjust the pressure field within a few model time steps.)

### b. Results

The simulated vortex quickly changes structure, especially near the surface, at the beginning of this simulation (hereafter referred to as the “baseline” simulation). A new quasi–steady state emerges after roughly 1000 s (Fig. 5). Most notably, the height of maximum wind speed is lower than in the axisymmetric simulation, the maximum value of wind speed is larger, and the periodic inertial oscillation above the boundary layer is no longer apparent (cf. Figs. 3 and 5).

Instantaneous velocity fields at are shown in Figs. 6 and 7. Vertical cross sections (bottom panels) show clearly turbulent flow in the inner core of the simulated vortex, as in previous studies of tornado-like vortices (e.g., Lewellen et al. 1997, 2000). Maximum horizontal and vertical wind speeds are roughly a factor of 2 larger than in the axisymmetric simulation near the surface ( m). These relatively large wind speeds are not necessarily attributable to localized gusts; flow near the surface can be nearly circularly symmetric (as shown in Figs. 6 and 7) and the near-surface structure is often quasi steady, although there are occasional oscillations in the structure and intensity for this case.

In contrast to the clearly turbulent character of the simulated flow in the inner core, there is little evidence of resolved turbulent fluctuations for radii greater than roughly 1 km, despite sufficient resolution out to a radius of 2 km. As discussed in the introduction, by formulation our LES model is supposed to produce turbulent fluctuations on the model grid (except, of course, near the surface, where is enforced, and eddies are typically too small to be resolved). Incorrect mean flow fields can develop in LES models that do not contain any turbulent fluctuations; for example, excessive vertical wind shear typically develops near the surface of neutral shear-driven atmospheric boundary layers (e.g., Markowski and Bryan 2016).

To analyze the properties of the resolved turbulent fluctuations in these simulations, we follow a common procedure: we define a mean state (denoted by angled brackets, e.g., for a generic variable *α*) and instantaneous fluctuations from that average (denoted by prime superscripts, e.g., ). However, because our simulated flows feature horizontally *heterogeneous* conditions associated with the spatially varying tornado-like vortex, some different techniques are required as compared to the more-common numerical studies of atmospheric boundary layers that have horizontally *homogeneous* conditions. In our case, an azimuthal average on a cylindrical grid might seem like an obvious choice. However, at small radii the number of points in each average can be quite small, and sometimes only one or two turbulent eddies are included in such averages. After evaluating several methods, we chose instead to define the mean state based on time averages, which is feasible here because the simulated vortex is roughly in statistically steady state, and is nearly stationary on the model grid. For the following analyses, we use 3D output every 10 s from s. The native model variables, notably the Cartesian components of velocity, are used for this averaging. Then, perturbations are calculated at every grid point, which can then be used to calculate products of perturbations (e.g., where *β* denotes a second generic variable). The products of perturbations are then averaged in time to determine resolved turbulent fluxes. For example, the resolved turbulent flux of in the *z* direction is , which is a three-dimensional diagnostic variable. For analysis purposes, these derived fields are then interpolated to a cylindrical grid, and are then averaged azimuthally.

Resolved turbulence kinetic energy (TKE) is defined as . The azimuthally averaged field of is shown in Fig. 8a. Maximum values, which exceed 1000 m^{2} s^{−2}, occur in the corner-flow region at small radius and height. In general, decreases as either *r* or *z* increases in the tornado core region of this simulation.

To aid the analysis of the inflow region (which, as shown in Figs. 6 and 7, seems not to produce turbulent fluctuations), vertical profiles of TKE are shown for two locations in Figs. 8b and 8c. The gray-shaded regions in these panels highlight layers in which the subgrid TKE () is greater than 10% of total TKE [i.e., ]. [Note that is a model-predicted variable, via (8)]. Generally speaking, the ratio is expected to be less than 0.1 for LES (i.e., most of the turbulent fluctuations are contained in resolved-scale motions); the exception, of course, is near the surface where turbulent eddies typically become very small (see section 2). The condition technically exists nearly everywhere in the fine-mesh subdomain of the baseline simulation, except below 25 m at m (Fig. 8b) and below 40 m AGL at m (Fig. 8c). These results are a bit surprising because turbulent fluctuations seem insignificant for km in the instantaneous fields shown in Figs. 6 and 7. We suspect that inertial waves, which are not the focus of this analysis, contribute to in these regions.

Most important to this study are turbulent fluctuations in the boundary layer of the tornado (roughly, the lowest few hundred meters AGL), which can transport high angular momentum downward leading to locally high wind speeds. Figure 9 shows an analysis of the magnitude of vertical turbulent flux of horizontal momentum, defined herein as

where the resolved component () is determined in the same manner as above (i.e., using time averages to define the mean state) and the parameterized/subgrid component () is given by (6). As shown in Fig. 9, is almost entirely parameterized in the boundary layer for m; that is, the resolved component is essentially not playing any role in the boundary layer inflow (which is of order 100 m deep, and is illustrated with a dashed line in Fig. 9a). Closer to the vortex core, the resolved component is larger in amplitude near the top of the inflow layer (yellow colors for m in Fig. 9a), but a large ratio of to (gray shading in Fig. 9c) extends up to 40 m AGL, or 16 grid points from the surface.

The very weak resolved turbulent fluctuations throughout the baseline simulation may be surprising considering the sufficient resolution and strong wind shear (in both radial and vertical directions) in the fine-mesh subdomain. However, this simulation features a transition from a completely parameterized turbulence scheme at large radii (2 km), and air flows inward toward the fine-mesh subdomain from this highly parameterized (and highly viscous) region.

Several numerical modeling studies of this type (i.e., LES embedded within a region of parameterized turbulence) have found that resolved turbulent fluctuations need to either be inserted or triggered continuously to facilitate the transition from completely parameterized to mostly resolved turbulence. Such techniques have a long history in computational physics. For example, the *fringe method* was developed by Spalart (1988). Related techniques were later developed by Lund et al. (1998), Mayor et al. (2002), and Inoue et al. (2014), among others. Some variant on these approaches could be used in the transition zone of our study. A barrier to these approaches for the present study, however, is that many of these methods implicitly assume that mean conditions are not too dissimilar, horizontally, across the domain. For example, in the “perturbation recycling” technique of Mayor et al. (2002) the turbulent fluctuations from a fully turbulent part of the domain are reinserted (i.e., recycled) into the upstream part of the domain. In our simulations, the mean wind speed and shear are considerably different near the center of our fine-mesh subdomain than at the transition zone, and thus turbulent fluctuations may not be appropriate for both locations.

Other techniques, such as a method that inserts random potential temperature perturbations in a transition region (e.g., Mirocha et al. 2014; Muñoz-Esparza et al. 2014), are designed to develop turbulence rapidly *within* the LES domain (without being “recycled” from another part of the domain). We considered using a version of such a technique, which is attractive for its simplicity and ease of implementation. However, we feel this approach can be problematic for flows with high wind speeds (e.g., tornadoes and hurricanes) because a large portion of the domain would be essentially devoted to developing turbulence. For example, we note that inflow velocity is *O*(5) m s^{−1} at our transition zone, and assuming an eddy turnover time of *O*(100) s, then a developing eddy will move *at least* 500 m into the fine-mesh subdomain before being fully developed (and, more likely, at least twice that distance). The computational burden of our application [including the rather small time step (of order 0.01 s) needed to maintain numerical stability for 100+ m s^{−1} flows] discourages us from devoting such a large portion of our computational resources for generation of turbulence. David Lewellen (2017, personal communication) has noted that the relatively small-scale turbulence in the surface layer has eddy turnover times s, and thus requires less distance from the transition zone to develop. The rough scale analysis above refers to a fully developed boundary layer, which extends above the surface layer.

There are other reasons we might be concerned about having a generation region for turbulence, such as how the “outer flow” vortex (above roughly 200 m) in our simulations might respond to changes in near-surface turbulence that would span the transition zone and generation regions. That is, the flow *above* the boundary layer, which varies radially and is part of the solution in our modeling studies (i.e., it is not a “given”), ideally should not reflect the essentially artificial changes in modeling strategy at the transition zone. For these reasons, we decided a rapid transition between parameterized and resolved turbulence in the transition zone was desirable for our study.

## 5. Precursor simulation 2: Doubly periodic LES

After testing several methods, we decided to use a method in which turbulent fluctuations develop in an offline “precursor” simulation that can be run for several thousand seconds, thereby ensuring a fully developed turbulent boundary layer. The velocity fields from this simulation are then inserted (i.e., “injected”) into the transition zone of the 3D vortex simulation. This type of technique has been used in many studies, primarily in the engineering literature (e.g., Thomas and Williams 1999; Keating et al. 2004), although several recent studies in the atmospheric sciences have utilized such a technique (e.g., Churchfield et al. 2012; Cunningham et al. 2013; Munters et al. 2016).

For convenience, our precursor LES is simply a “traditional” simulation of shear-driven turbulence on a fairly small Cartesian domain using doubly periodic lateral boundary conditions, similar to simulations by Moeng and Sullivan (1994) and Conzemius and Fedorovich (2006). An important difference from these previous studies, however, is the addition of tendency terms to the horizontal velocity equations to account for processes associated with rotating flow. Specifically, we add terms accounting for centrifugal acceleration, radial pressure gradient, and radial advection in the presence of a mesoscale horizontal gradient of velocity. The technique was developed and evaluated for simulations of tropical cyclone boundary layers by Bryan et al. (2017) and simply adds terms to the horizontal velocity equations of (5a):

For ease of interpretation, we assume the domain for precursor 2 simulations is located to the east of the tornado-like vortex; thus, the field is analogous to radial velocity, and the field is analogous to tangential velocity. Four parameters on the right-hand side of (13) are specified at the beginning of the simulation and are held fixed: the rotation rate , a reference radius *R*, a tangential velocity *V* (technically, a gradient wind), and the radial gradient of *V* (i.e., ). As noted by Bryan et al. (2017), the variables *V* and may be specified as vertical profiles (i.e., functions of *z*), but herein we set these parameters to constants, for simplicity. The remaining parameters on the right-hand side of (13) are calculated from the model-predicted fields, where tildes denote horizontal averages over the entire domain; for example, is the horizontal average of the component of velocity. As explained by Bryan et al. (2017), these “mesoscale tendency” terms account for the most important dynamical effects in a rotating flow. These terms are applied uniformly on each level of a Cartesian domain. Consequently, the simulated flow is not “curved” (as in the 3D baseline simulation). Nevertheless, the simulated mean wind profiles are remarkably similar to those from strongly curved, rotating flow.

For this “precursor 2” simulation we use a 2 km 4 km 1 km domain using the same grid spacing as before (5 m, 2.5 m). The governing equations and numerical methods are exactly the same as before except for the inclusion of the mesoscale tendency terms, in (13), and we set in (5a). The lower boundary condition is also the same as before (see section 2). The lateral boundary conditions are periodic in both directions, and a Rayleigh damper is included above 750 m AGL that relaxes the predicted velocity fields back to their initial values.

The initial conditions are and , plus random small-amplitude (0.25 m s^{−1}) perturbations to and . For this simulation we set km (the radius of the transition zone in the baseline simulation) and we use the same value of from that simulation. The values of *V* and come from the precursor 1 simulation at km; specifically, *V* = 10.4 m s^{−1} and = s^{−1} (which were obtained by averaging in the vertical from 0.5 to 2.0 km AGL).

Figure 10 shows instantaneous views of vertical velocity after 4000 s of integration. The simulated flow requires 1000 s to become fully turbulent and settle into a statistically steady state. The domain-mean fields (solid lines in Fig. 11a) feature essentially radial inflow (i.e., layer with ) in the lowest400 mwith a maximum amplitude of *O*(5) m s^{−1}, and a tangential velocity component () that is strongly sheared in the lowest 200 m; both of these features are qualitatively similar to mean fields in the precursor 1 simulation at km (dashed lines in Fig. 11a). The velocity fields from the two precursor simulations are not identical; for example, the “inflow” layer is deeper in the precursor 2 simulation. However, differences are expected for several reasons, including the following: 1) the updraft forcing term *F* is set to zero in the precursor 2 simulation, but is influencing the flow in the precursor 1 simulation (and is responsible for the m s^{−1} values of for m); 2) a tunable length scale is used for parameterized turbulence in the precursor 1 simulation, and the value chosen herein (5 m) is probably too small; and 3) large-scale subsidence has been neglected in the precursor 2 simulation, for simplicity, but this effect exists in precursor 1 simulations and acts to limits the growth of the boundary layer. Nevertheless, the two wind profiles have the same qualitative structures, and are quite similar in the surface layer (roughly, the lowest 50 m AGL in this case). Also, most important for the present purpose is that the domain-average parameterized and resolved TKE in the precursor 2 simulation (Fig. 11b) demonstrate the desired properties for our study [i.e., most of the TKE is contained in the resolved scales (except near the surface)].

During this simulation, we write perturbation velocity fields to a file every other time step along a 2D (*y*–*z*) slice in the middle of the domain (i.e., at km). Specifically, we write , , and on a 4 km 0.5 km slice (i.e., 2D cross section); we only write the lowest half of the domain (0–0.5 km AGL) to reduce file size, and because fluctuations above 0.5 km are relatively weak. These fields are then used to inject turbulence into our 3D domain, as described below.

## 6. A 3D simulation with eddy injection

### a. Formulation

As explained in section 1, a primary goal of our simulations is to have fully developed turbulent perturbations in the boundary layer inflow of the tornado-like vortex. The simulation developed in this section, utilizing the technique we refer to as eddy injection, is identical to that in section 4 except for a set of tendency terms that are only applied in a narrow ring within the transition zone (gray region of Fig. 4). Specifically, we take the perturbation fields , , and written from the precursor 2 simulation (see section 5) and “nudge” them into the 3D domain in this ring (as described below). These velocities are assumed to be from a cylindrical coordinate system, with being the radial component and being the tangential component. Because the domain from the precursor 2 simulation is only 4 km long, but the transition zone is km long, we simply replicate the precursor 2 results, utilizing its periodic lateral boundary condition, as needed to circle the entire transition zone. Although there is technically a discontinuity in conditions where this procedure begins/ends (because the transition-zone length of 12.6 km is not exactly divisible by the precursor 2 domain length of 4 km), this situation does not become problematic because we nudge the velocity fields toward conditions from the precursor 2 simulation (i.e., we do not force the exact conditions into the transition zone).

Specifically, we apply a relaxation term that nudges the actual perturbation velocities at each grid point (relative to the azimuthal average) to the value from the precursor 2 simulation:

where, as before, subscripts *r* and *ϕ* denote radial and tangential components of velocity, respectively; and curled brackets denote azimuthal averages. The *γ* term varies between 0 and 1, and is included to turn off the perturbation tendencies above a certain height. (As explained above, the perturbation file is only written from the surface to 500 m AGL.) Here, we use , with m and m. The relaxation time scale is set to 1 s. The cylindrical-coordinate velocity tendencies in (14a) and (14b) are converted into Cartesian components and are added to the CM1 Cartesian velocity equations at grid points in the aforementioned ring, which is located in the middle of the transition zone (see Fig. 4). The ring over which these perturbations are applied is only two grid points wide; specifically, these tendencies are only applied at grid points that are located within of the center of the transition zone.

This procedure is conducted every time step. Because the time steps of the precursor 2 simulation are typically different from the 3D LES (which uses an adaptive time step described in section 2) we simply use the closest relevant time from the precursor 2 simulation for the perturbation fields (the double-primed variables above); interpolation between neighboring times is likely unnecessary because fluctuations on the scale of a few time steps are negligibly small in CM1 output (e.g., Worsnop et al. 2016, manuscript submitted to *Bound.-Layer Meteor*.).

### b. Results

A snapshot of horizontal velocity is shown in Fig. 12. Compared to the baseline simulation (Fig. 6) the character of the simulated flow is quite different. Of particular note is the development of “streaks” of relatively high and low wind speed near the surface (Fig. 12a). As expected, these streaks align with strips of rising and sinking motion associated with quasi-2D roll circulations (Fig. 13a). Similar features were noted in simulations by Lewellen and Zimmerman (2008). Of particular importance, these resolved turbulent eddies near the surface extend all the way to the transition zone in this study, denoted by the dashed lines in Figs. 12 and 13.

Average velocity fields are shown in Figs. 14a–c along with differences from the baseline simulation in Figs. 14d–f. Radial velocity is primarily weaker and deeper with eddy injection. The mean azimuthal velocity profile is less sheared below 100 m AGL with eddy injection, which can be seen in Fig. 14e as relatively larger values of below 25 m with weaker values from roughly 25 to 100 m. The mean vertical velocity field shows the greatest structural difference between the two simulations: maximum is nearly on the central axis () in the baseline simulation (Fig. 5c), but is clearly located off the central axis with eddy injection. The central downdraft is stronger, and reaches closer to the surface, with eddy injection.

These differences in mean flow are not a direct consequence of the eddy injection method, which inserts perturbations only; thus, on average over time, the method adds no net tendency to the three velocity components. Rather, the different mean fields develop *within* the LES subdomain as a consequence of the different representations of turbulence. Figure 15 illustrates how the differences develop. Average tangential velocity is shown at four locations in Figs. 15a–d. On either side of the transition zone, is essentially identical in the two simulations (Figs. 15c,d). Inward of the transition zone (Figs. 15a,b), the baseline profiles (black curves) become more sheared than in the simulation with eddy injection (red curves).

Figures 15e–h show average profiles of the vertical turbulent stress (including both resolved and parameterized components) for the azimuthal direction (), which explain the primary reason for the differences in . Outside the transition zone, the two simulations have essentially identical profiles (Fig. 15h). Just inward of the transition zone, the profiles are substantially different (Fig. 15g), with the stress in the simulation with eddy injection being substantially more negative in the lowest 100 m AGL. As shown below, the primary difference is attributable to the lack of resolved turbulent eddies in the baseline case. The profiles remain different at regions farther inward (Figs. 15e,f) and are consistent with the relative changes in (noting that ).

The magnitudes of and [see (12)] for this simulation are shown in Fig. 16. Values of are much larger in most of the domain compared to the baseline simulation, especially outside the core of the vortex (roughly, m) (cf. Figs. 9a and 16a). The parameterized component is less than 10% of the total value above roughly 20 m AGL in the simulation with eddy injection (Figs. 16c,d). Furthermore, is more than an order of magnitude larger in the simulation with eddy injection over the lowest 100 m AGL. Most important, the radial distribution of (=) is relatively smooth across the transition zone in the simulation with eddy injection, but abruptly changes (in magnitude and extent from the surface) in the baseline simulation due to the near absence of a resolved component in that case (Figs. 9a,b). Readers interested in further details of the momentum budget are referred to section 4 of Nolan et al. (2017).

Analysis of TKE is shown in Fig. 17. With eddy injection, the contours of resolved TKE () extend all the way to the edge of the fine-mesh subdomain at km; in contrast, is much smaller in the baseline simulation for (cf. Figs. 8a and 17a). The difference field (Fig. 17d) clearly illustrates these differences.

Comparison of profiles of and show fairly similar results at m (cf. Figs. 8b and 17b), although the layer with (gray shading) is roughly half as deep in the simulation with eddy injection. Differences at m (cf. Figs. 8c and 17c) are more pronounced; near m (dotted line) the resolved TKE is an order of magnitude larger with eddy injection.

As another example of the differences in these simulations, Fig. 18 shows an analysis of the maximum instantaneous wind speed over the final 1000 s of these simulations. Specifically, for every Cartesian grid point in the CM1 domain we search for the maximum instantaneous value of horizontal wind speed (using 3D output files between and s, which are written every 10 s); then, we interpolate the results to a cylindrical grid and search for the maximum value at every radius and height. The result (Fig. 18) shows that maximum instantaneous gusts in the baseline simulation are substantially stronger ( m s^{−1}) than in the simulation with eddy injection. Near the location of strongest time-averaged wind speed [(*r*, *z*) (100, 50 m); see white crisscross () symbol] the maximum gust is more than 100 m s^{−1} weaker in the simulation with eddy injection (Fig. 18c). Similar results are found for maximum *vertical* velocity (Fig. 19). If debris were included in these simulations, we suspect there would be a substantial difference in the distribution of lofted debris, considering these large differences in near-surface vertical velocity.

We have repeated the analyses from this section for simulations with other input parameters (e.g., vertical velocity forcing, rotation rate, and surface roughness length) and have reached the same overall conclusions. As an example, the instantaneous horizontal wind speed for a smaller swirl ratio (, a factor of 4 smaller than before) is shown in Fig. 20. In this case, the simulated vortices have a “single-celled below and double-celled above” structure (e.g., Rotunno 2013) typical of moderate swirl ratio, although vortex breakdown occurs at a lower altitude for the simulation with eddy injection. Most important, the simulated flow in the boundary layer of the baseline configuration has very weak resolved turbulent fluctuations (Fig. 20a). When eddy injection is included, the streaks associated with quasi-two-dimensional rolls are apparent (Fig. 20c), and the boundary layer inflow is deeper (not shown), similar to the relatively higher swirl vortices analyzed earlier.

Our overall conclusion from such comparisons is that the simulated vortices are substantially different when eddy injection is included. However, drawing a simple set of conclusions about *how* the vortices are different turns out to be tricky. When eddy injection is included, the flow at some distance from the tornado core (roughly, for m) seems to be more like low-swirl and/or low–Reynolds number cases from our axisymmetric model study (RBND), in the sense that boundary layer inflow becomes deeper and weaker. However, in the inner core, when eddy injection is included the simulated vortices generally behave like higher-swirl and/or higher–Reynolds number cases from that study, in the sense that vortex breakdown happens at a lower height (e.g., Fig. 20; also, cf. Figs. 5 and 14). Further analysis of this topic, using angular momentum budgets, is included in a companion paper (Nolan et al. 2017), which supports the latter interpretation of higher “effective” swirl ratio in the inner core for simulations with eddy injection.

The maximum value of wind gusts is also somewhat complex. Overall, we find that peak horizontal and vertical wind speeds tend to be lower when eddy injection is included, as shown in Figs. 18 and 19. However, relatively low- cases suggest the opposite conclusion, at least near the surface (cf. Figs. 20a and 20c). A partial explanation for this complexity depends on which simulation is closest to the “optimal vortex” with an end-wall vortex structure (e.g., Fiedler and Rotunno 1986; RBND).

## 7. Summary

A framework to study wind speeds in tornado-like vortices has been described in this article. Following a companion article (RBND), the “Fiedler Chamber” type of approach is used in which updraft forcing is specified in the middle of a closed, rotating domain. A fine-mesh subdomain is used in the center of the overall domain to conduct large-eddy simulation (LES). A simple but effective parameterized turbulence scheme is used for the rest of the domain. To reduce computational cost, an axisymmetric version of the model is first integrated for *O*(10 000) s to spin up the flow to an approximately steady state. In a baseline 3D simulation, these initial conditions, plus random perturbations, are then integrated for several thousand seconds.

The results of this baseline simulation are mixed, in our view. The inner core of the tornado-like vortex clearly produces resolved turbulent fluctuations (as expected for the design of these simulations); however, the boundary layer inflow does not clearly contain turbulent fluctuations. Simulated flows from LES models that do not have resolved turbulent fluctuations can have erroneous characteristics, such as excessive vertical wind shear (e.g., Markowski and Bryan 2016). Further analysis of the baseline simulation shows that large percentages of both the turbulence kinetic energy (TKE) and the vertical turbulent flux of horizontal momentum are associated with the subgrid turbulence parameterization. We suspect that some recent LES studies of tornado-like vortices (e.g., Kuai et al. 2008; Maruyama 2011; Natarajan and Hangan 2012; Liu and Ishihara 2015; Bodine et al. 2016) have the same characteristics, based on figures within those articles.

Several methods have previously been developed to address the underlying problem, which is the lack of development (or slow development) of resolved turbulent fluctuations near the interface between flow with completely parameterized turbulence and flow that has resolved turbulent fluctuations (e.g., Lund et al. 1998; Mayor et al. 2002; Keating et al. 2004; Muñoz-Esparza et al. 2014; Munters et al. 2016). We decided to develop turbulent fluctuations in a separate “precursor” LES, which uses periodic lateral boundary conditions on a smaller domain, using the technique of Bryan et al. (2017) to account for a rotating environment. This precursor simulation is integrated for sufficient time (1000 s) to develop a turbulent boundary layer that is in a statistically steady state. Perturbation wind fields from a two-dimensional slice of this simulation are then injected into a 3D simulation in a ring near the edge of the fine-mesh subdomain. The boundary layer clearly contains turbulent fluctuations in the subsequent simulation, both near the injection ring, and throughout the boundary layer inward of the ring. Several aspects of the simulated vortex are different, compared to the vortex in the baseline simulation; most notably, the subgrid turbulence scheme is responsible for a small percentage of the vertical turbulent momentum flux and TKE (except, as expected, very close to the surface), the boundary layer inflow is deeper, and several measures of wind speed usually have lower values in the simulation with eddy injection.

Further analyses of simulated tornado-like vortices using this modeling framework are presented in two companion papers (Nolan et al. 2017; Dahl et al. 2017). Future work on the framework described herein could look into optimized mesh sizes: the relatively large (4 km 4 km 1 km) fine-mesh subdomain used herein is probably larger than necessary for most applications. We chose this subdomain size somewhat arbitrarily, albeit with the goal of ensuring that turbulence in the boundary layer inflow has sufficient time to mature before reaching the inner-core region. Also, the depth of the fine-mesh subdomain could probably be reduced for most cases. However, it is possible that higher resolution in the surface layer, such as grid spacing m (horizontal and vertical) in the lowest 10 m AGL, might be needed to properly resolve turbulent features in some tornadoes. High-resolution observational data, with similar resolution, would also be beneficial.

Sensitivity experiments should also be conducted with different LES subgrid schemes, and particularly schemes that are not entirely dissipative (i.e., the effects of backscatter from subgrid scales should be studied). Subgrid models that better handle the transition from parameterized to resolved turbulence near the surface (e.g., Mason and Thompson 1992; Sullivan et al. 1994; Kosović 1997; Chow et al. 2005; Brasseur and Wei 2010) should also be evaluated.

Finally, this initial study focused intentionally on the simplest case (i.e., a quasi-steady and stationary vortex without asymmetries in the environment or underlying surface). In a companion paper (Dahl et al. 2017), we used the eddy injection technique for simulations of translating vortices, which develop notable asymmetries. We use only one precursor 2 simulation for each simulation, although it is possible that the mean conditions might become substantially different along the ring where the eddies are injected. Thus, perhaps different eddy fields, from different precursor 2 simulations, could be injected in different regions around the vortex, relative to the translation direction. It is unclear, at this time, whether any further modifications to the technique described herein would be needed in this case. Future studies might consider such a modification of the technique for other applications, such as simulations with heterogeneous surface conditions (e.g., land use and/or terrain).

## Acknowledgments

D. Nolan and N. Dahl were supported by the National Science Foundation through Grant AGS-1265899. The authors acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc), project UMIA0009, provided by NCAR’s Computational and Information Systems Laboratory. The authors appreciate the helpful comments and suggestions from Dr. David Lewellen and an anonymous reviewer.

## REFERENCES

*An Introduction to Fluid Dynamics.*Cambridge University Press, 615 pp.

*24th Conf. on Severe Local Storms*, Savannah, GA, Amer. Meteor. Soc., 8B.1. [Available online at https://ams.confex.com/ams/pdfpapers/141749.pdf.]

*Encyclopedia of Atmospheric Sciences*, 2nd ed. G. North, F. Zhang, and J. Pyle, Eds., Academic Press, 232–240.

*Turbulence in the Atmosphere.*Cambridge University Press, 393 pp.

## Footnotes

^{a}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).