Introduction

The effects of surface drag are important throughout the atmospheric boundary layer (ABL). Surface drag, in addition to baroclinicity, plays a role in the development of the highly sheared ABLs that are favorable for the formation of severe convective storms. Moreover, surface drag has been found to influence the dynamics of convective storms themselves. For example, Schenkman et al. (2012) found that including surface drag in a numerically simulated mesoscale convective system (MCS) promoted the generation of a horizontal rotor that ultimately abetted the development of an intense vertical vortex. In a bow echo simulation by Xu et al. (2015), it also was concluded that surface drag was an important source of circulation for the eventual vertical vortices, although the contribution of surface drag was obtained as a residual in the circulation analysis. In simulations of supercell thunderstorms, horizontal vorticity generation by surface drag, with subsequent tilting and stretching, has been implicated in the development of some tornado-like vortices (Schenkman et al. 2014;Roberts et al. 2016;Yokota et al. 2016). Furthermore, in addition to surface drag being a source of horizontal vorticity, it is also well known that surface drag can intensify a vertical vortex by enhancing radial inflow within the boundary layer of the vortex, thereby enhancing the convergence of angular momentum (e.g., Burggraf et al. 1971;Lewellen 1976Lewellen , 1993;;Davies-Jones et al. 2001;Rotunno 2013). However, the idealized vortex simulations of Nolan et al. (2017) and the idealized supercell simulations of Roberts et al. (2020) indicate that a tornado's strength may increase with surface roughness only to a certain point.

In atmospheric numerical models, regardless of whether they use Reynolds-averaged Navier-Stokes (RANS) or large-eddy simulation (LES) approaches, the surface drag is computed by models of near-surface turbulence (also known as surface-layer parameterization schemes or wall models) that use wind vectors at the lowest one or more grid levels as inputs. The aforementioned convective-storm flows are frequently characterized by substantial unsteadiness and horizontal heterogeneity. Convective storm modeling studies, however, have employed models of near-surface turbulence derived for statistically steady and horizontally homogeneous conditions. Such approaches are known as equilibrium approaches.

The Monin-Obukhov similarity theory (MOST; Monin and Obukhov 1954) has been the most commonly used equilibrium approach. MOST assumes an instant, local relationship between near-surface horizontal velocity components and the surface shear stress; this is known as the "semi-slip" lower boundary condition. Note that MOST is an extension of the law of the wall (LOTW; Prandtl 1933) from neutral to nonneutral conditions, whereas the LOTW is derived for statistically steady, horizontally homogeneous, channel flows. Unsurprisingly, large deviations from MOST have been observed within convective storm outflow and in near-storm environments (Markowski et al. 2019). A recent LES study shows that accounting for the evolution of near-surface turbulence along curved trajectories (which is no longer an equilibrium approach) may intensify tornadoes, particularly during time periods in which tornadoes are highly unsteady (Wang et al. 2020). These observational data and LES results raise concerns about simulating convective storms with near-surface turbulence modeled using equilibrium approaches.

The primary objective of this work is to investigate the potential consequences of the assumption of statistically steady and horizontally homogeneous near-surface conditions, which has been used in convective storm simulations via the semi-slip lower boundary condition. In addition to the nonequilibrium approach proposed by Wang et al. (2020), an engineering nonequilibrium approach known as the thin-boundary layer equationfoot_0 (TBLE) is tested. The TBLE was first proposed by Balaras et al. (1996) for LES of smooth-wall boundary layers where roughness elements are much smaller than the viscous length scale. The unsteadiness and horizontal heterogeneities of near-surface turbulence are accounted for by solving RANS equations simplified for a "wall layer" between the surface and the lowest LES grid level for horizontal velocity components. Balaras et al. (1996) used TBLE to improve the reproduction of flows in the corner regions of a square duct. Cabot (1996) refined the eddy-viscosity model employed by TBLE and further improved the reproduction of flow separation downstream of a step. These successful case studies make TBLE an attractive candidate for convectivestorm applications.

In this work, the TBLE framework is implemented into a widely used atmospheric model, the Cloud Model 1 (CM1; briefly described in section 2a). The implementation details, including modifications needed to adopt TBLE for atmospheric applications over roughness elements much larger than the viscous length scale (i.e., fundamentally different from smoothwall boundary layers), are described in section 2b. In sections 3 and 4, CM1 LES runs are performed for an idealized ABL and an idealized tornado, respectively, using four different models of the near-surface turbulence: (i) the semi-slip lower boundary condition (which is an equilibrium approach); (ii) the nonequilibrium approached proposed by Wang et al. (2020) (hereafter referred to as the turbulence-memory model); (iii) a TBLE with the eddy-viscosity model used by Balaras et al. (1996), which is referred to as Balaras TBLE hereafter; and (iv) a TBLE with the eddy-viscosity model used by Cabot (1996), which is referred to as Cabot TBLE hereafter. An idealized ABL is simulated for code verification purposes, while an idealized tornado is simulated to study the sensitivity of a rapidly evolving dynamic system with strong horizontal heterogeneities to the representation of near-surface turbulence. Conclusions are presented in section 5.

Methodology

a. CM1 LES

The LES are performed using CM1 (version 19.6), a compressible and nonhydrostatic model. For simulations of a dry atmosphere, CM1 solves the prognostic equations of velocity, potential temperature, and nondimensional pressure (see the appendix of Bryan and Morrison 2012). The time integration uses a third-order Runge-Kutta scheme. A fifth-order upwind advection scheme is used, and a fifth-order weighted essentially nonoscillatory (WENO) scheme is applied to advection of scalars (as recommended by Wang et al. 2021) on the final Runge-Kutta step. The subgrid-scale (SGS) momentum and scalar fluxes are computed using a turbulent kinetic energy (TKE) closure (Deardorff 1980), which solves an additional prognostic equation for the SGS TKE.

The lower boundary conditions for horizontal momentum equations are computed using various models of near-surface turbulence. The semi-slip lower boundary condition is already an available option in CM1, the turbulence-memory model has been recently proposed to account for the evolution of turbulence along curved air parcel trajectories (Wang et al. 2020), and the TBLEs are new models implemented into CM1 (explained in section 2b).

b. Implementing TBLEs into CM1

A TBLE takes as inputs the resolved fields computed within an LES and calculates unresolved fluxes by solving RANS equations within a wall layer between the surface and the lowest LES grid level for horizontal velocity components. For a dry atmosphere, the full three-dimensional RANS equations, employing the Boussinesq approximation in a rotating frame of reference, are given by (Holton and Hakim 2013, chapter 8.1.1):

(1)

Here an overbar represents an ensemble average, and a superscript prime represents a perturbation from the ensemble averaged value; u, y, and w are zonal, meridional, and vertical velocities; u 0 and r 0 are hydrostatic-base-state potential temperature and density, respectively; and u and p are perturbations of potential temperature and pressure, respectively, from their hydrostatic base states (meaning that u and p are their ensemble-averaged "perturbation from hydrostatic base state"); f is the Coriolis parameter; n is the molecular kinematic viscosity; and n h is the molecular thermal diffusivity.

When the ensemble-averaged flow field is predominantly horizontal, the horizontal scales are much larger than the vertical scales (according to the continuity equation), meaning that turbulence statistics vary mainly in the vertical direction. This assumption is applicable to most near-surface flows in the atmosphere except for some extreme cases (e.g., near the core of a tornado). Applying this assumption to (1) yields

(

If buoyancy is negligible in the w budget, and the vertical turbulent transport terms are leading-order terms in both u and y budgets, then scale analysis suggests a main balance for the w budget (Pope 2000, chapter 5.2.1):

Taking the horizontal derivative of (3) yields

where the spatial variation of r 0 (which is nonzero only in the vertical direction) is negligible compared to the spatial variations of other variables, according to the Boussinesq approximation. Integrating (4) from the surface to an arbitrary height (z) within the wall layer yields

where p sfc is the perturbation pressure at the surface. Thus, when the horizontal variations of w w are negligible compared to the horizontal perturbation pressure gradient force (HPPGF), the HPPGF remains approximately constant with height. The vertically constant HPPGF conforms with earlier boundary layer work (e.g., Batchelor 2000, p. 305) and is also assumed by equilibrium approaches (Pope 2000, chapter 7.2.1). This key assumption taken by the TBLE needs to be examined for each of the cases of application (see sections 3b and 4b).

Combining (1) through (5) with the above assumptions of horizontally homogeneous turbulence and neutral stratification near the surface (which renders the internal energy equation irrelevant), the RANS equations within the wall layer can be simplified as (cf. Balaras et al. 1996, section 2):

Eddy-viscosity models are used for the turbulent fluxes:

where the eddy viscosity (n t ) has a dimension of a length scale multiplied with a velocity scale. At least two models of n t have been used in previous TBLE applications. Balaras et al. (1996) used a mixing-length-type eddy-viscosity modelfoot_4 :

where k 5 0.4 is the von Kármán constant, z 0 is roughness length, k(z 1 z 0 ) provides the characteristic length scale, and k(z 1 z 0 )|S| provides the characteristic velocity scale. Here |S| 5 (2S ij S ij ) 1/2 is the magnitude of the mean strain-rate tensor, where

Given the wall-layer assumption that velocity statistics vary much more rapidly in the vertical direction than in the horizontal directions, the magnitude of the mean strain-rate tensor can be approximated as

Cabot (1996) used a different approach to model the eddy viscosity:

where the friction velocity u * provides the characteristic velocity scale, and k(z 1 z 0 ) provides the characteristic length scale. The major difference between (8) and ( 11) is whether the characteristic velocity scale is determined by the local mean shear (Balaras et al. 1996) or the wall stress (Cabot 1996). Although both eddy viscosities given by ( 8) and ( 11) are able to recover LOTW within the inertial sublayer embedded within a statistically steady, fully developed, neutral boundary layer above a horizontally homogeneous surface, using the local mean shear to determine the characteristic velocity scale as suggested by ( 8) becomes unphysical when a velocity profile involves a local extremum within the wall layer.

When being used as an approach of modeling near-surface turbulence in an LES, the HPPGF terms in ( 6) take the values computed by the LES at its lowest grid level for horizontal velocity components (Balaras et al. 1996;Cabot and Moin 2000;Wang and Moin 2002). The top boundary conditions for u and y in (6) also come from values at their lowest LES grid level. The TBLEs use a horizontal grid spacing identical to that of the LES and a vertical grid spacing much smaller than that of the LES (illustrated by blue lines in Fig. 1).

All previous TBLE applications (to the authors' knowledge) have focused on smooth-wall boundary layers, where the roughness elements are much smaller than the viscous length scale. These TBLE applications have used a vertical grid spacing on the order of the viscous length scale, and the surface friction is given by the molecular viscous stress within the laminar viscous sublayer. Within a rough-wall boundary layer like the ABL, the roughness elements are much larger than the viscous length scale, and the surface friction is primarily contributed by pressure drag exerted by roughness elements. Because molecular viscosity is no longer important in a rough-wall boundary layer, using a TBLE vertical grid spacing as small as the viscous length scale is physically meaningless. Theoretically, the minimum TBLE vertical grid spacing is on the order of the roughness length (z 0 ), and the surface friction needs be computed as the turbulent shear stress above the roughness elements. As a preliminary attempt, the surface shear stress is computed using the LOTW with inputs taken from u and y values at a TBLE grid level which satisfies z . . z 0 (a condition required for the LOTW to be applicable, see Pope 2000, chapter 7.2.2).

Finally, because molecular viscosity is unimportant in a rough-wall boundary layer, the TBLE coupled with CM1 LES solves a further simplified set of RANS equations:

where the viscous forces in the u and y budgets in ( 6) have been neglected.

Testing the TBLE implementation by simulating an idealized atmospheric boundary layer

Theoretically, the mean-wind profile in a statistically steady, horizontally homogeneous, fully developed, neutral surface layer follows LOTW, which provides a means of validating LES code. Here simulations of a quasi-steady, horizontally homogeneous, fully developed, inertial-driven neutral ABL are performed to test: (i) whether the TBLEs are coded into CM1 LES correctly and (ii) whether the new application of TBLEs to a rough-wall boundary layer is appropriate.

a. Simulation configuration

The simulations employ a 10 km 3 10 km 3 4 km domain with periodic lateral boundary conditions, a free-slip upper boundary condition, and lower boundary conditions computed using models of the near-surface turbulence, including the semislip boundary condition, the recently proposed turbulencememory model (Wang et al. 2020, applied here with the model parameter g set to 0.1), and the TBLE approach with two different eddy-viscosity models described in section 2b. The LES grid spacing is 25 m in each direction, and the roughness length is z 0 5 0.12 m. When using TBLEs, the wall layer (i.e., z # Dz/2, where Dz is the LES grid spacing) is discretized into 20 layers, meaning that the TBLE vertical grid spacing is aboutfoot_5 0.625 m. The u and y values at the fifth TBLE grid level (i.e., z 5 z wl ≃ 3.125 m . . z 0 ) are taken as inputs to the LOTW to compute the surface shear stress.

Each simulation starts with a neutrally stratified ABL where u 0 5 300 K below z 5 1 km and increases at a rate of 1 K km 21 above z 5 1 km, with a random perturbation within 60.25 K added at each grid point to trigger turbulence. Rayleigh damping is applied above z 5 2.4 km. The initial velocity field is ũ 5 10 m s 21 and ỹ 5 w 5 0 everywhere, where the tilde represents spatial filtering over an LES grid. No Coriolis force is applied. The horizontal mean pressure gradient is zero, and therefore the flow is driven by inertia. Note that an inertia-driven ABL yields slowly decaying turbulent flow, which may not be ideal to reproducing LOTW derived for strictly statistically steady conditions. Nevertheless, the scope of the current section is to examine the TBLE implementation into CM1 LES, rather than reproduce LOTW. Here the simulation using a semi-slip lower boundary condition is exactly the same as the "SWENO-DX25DZ25" run in Wang et al. (2021), except for using a different version of CM1. Each simulation runs for 6 h with a time step of 0.4 s. Each simulation reaches a quasi-steady state at t 5 3 h or earlier, when the turbulence statistics normalized by the friction velocity no longer vary with time. The Reynolds-averaged results are estimated by taking the horizontal mean and then averaged over t 5 5-6 h.

b. Examining the TBLE assumption about HPPGF

The simulation with a semi-slip boundary condition is used to examine the assumption of a vertically constant HPPGF within the wall layer in order to understand whether this assumption is valid when the TBLE approach is used. Assuming vertically constant HPPGF is based on the assumption that the leading terms in the w budget are the vertical perturbation pressure gradient force (VPPGF) and vertical turbulent FIG. 1. Illustration of the wall layer for a TBLE within the LES domain and how surface shear stress is derived. Blue lines and fonts represents the TBLE part. Here z wl is the reference height in the wall layer where u and y are taken as input to the LOTW, the tilde over variables represents spatial filtering over an LES grid cell, and the overbar indicates the RANS variables in the wall layer. mixing (VTURB), as shown by (3) through (5). To investigate the validity of the assumed w budget near the surface, the w budget at z 5 Dz (the lowest w level above the surface) is examined. Figures 2a-f reveals that the leading terms in the w budget are not VPPGF or VTURB expected by (3), but the tendency (WTEND) and horizontal advection (HADV) terms. Unlike the RANS equations where turbulent motions are included in the turbulent flux terms, LES resolves turbulent motions, and HADV and WTEND are of larger magnitudes compared to the VPPGF. This violates the assumption given in (3), which in turn invalidates the assumption of vertically constant HPPGF in the wall layer. As a result, using the HPPGF calculated at the lowest LES grid level of horizontal velocities (z 5 Dz/2) for the entire wall layer is questionable.

To confirm the invalidity of assuming vertically constant HPPGF in the wall layer, the difference between HPPGFs at the lowest LES grid level and at the surface is estimated. Because the VPPGF is not a leading-order term in the w budget, one cannot solve for the surface perturbation pressure (p sfc ) by analytically integrating the w budget equation from the surface to the lowest LES grid level. In other words, p sfc can only be estimated by vertically integrating the VPPGF computed by the LES, which is equivalent to some mathematical extrapolation. By default, the CM1 code outputs psfc computed using a third-order extrapolation: FIG. 2. Histograms of terms in the w budget at z 5 Dz (estimated using the w budget in the LES domain) in the LES of an idealized ABL, including (a) vertical perturbation pressure gradient force (VPPGF), (b) w tendency (WTEND), (c) horizontal turbulent flux (HTURB), (d) vertical turbulent flux (VTURB), (e) horizontal advection (HADV), and (f) vertical advection (VADV). Here WTEND is obtained as the summation of all other terms. (g) The distribution of the difference of the HPPGFs at z 5 Dz/2 and at surface in (15) scaled by the HPPGF at z 5 Dz/2 in ( 14). The probability distributions are computed using results at t 5 6 h. psfc 5 15 p0:5Dz 2 10 p1:5Dz

where the subscripts indicate the heights of associated LES variables. The magnitude of HPPGF at the lowest LES grid level is given by

and the difference between HPPGFs at the surface and at the lowest LES grid level can be estimated as

where the vertical variation of density is neglected. Figure 2g shows that the mean DHPPGF/HPPGF 1 is 0.5, so DHPPGF is nonnegligible, violating the assumption of vertically constant HPPGF used in the TBLE. The question remains as to how much of an impact the violation of this assumption may have on the TBLE simulation results.

c. Results

Results of surface-layer mean shear are evaluated against LOTW filtered over the LES vertical grid spacing (the gray dotted line in Fig. 3; see more details in Wang et al. 2021). Specifically, the filtered LOTW is calculated as

and the filtered nondimensional shear is given by

Figure 3 shows that all LES runs overpredict the nondimensional shear at the second and the third LES grid levels, where the SGS model performance is challenged by underresolved energy-containing turbulent motions (e.g., Mason and Thomson 1992;Sullivan et al. 1994;Porté-Agel et al. 2000;Bou-Zeid et al. 2005;Yang et al. 2017;Wang et al. 2021). Compared to the semi-slip boundary condition and the turbulence-memory model, the TBLEs slightly amplifies the overprediction of the nondimensional shear (less than 10% difference). The slightly worse performance of the TBLEs is caused by the direct match of LES variables and RANS variables at z 5 Dz/2, the lowest level of the LES grid and the uppermost level of the wall layer for horizontal velocities. RANS models all scales of turbulence, whereas LES resolves the energy-containing turbulence, so using LES variables for RANS equations leads to overestimation of the momentum and the resulting surface shear stress (Wang and Moin 2002;Park and Moin 2014). Because TBLEs are used later to explore the influence of inhomogeneity, their primary value is not necessarily to better reproduce LOTW.

To examine the potential sensitivity of TBLE results to the variation of HPPGF with height within the wall layer, an additional Cabot TBLE simulation is performed with an HPPGF linearly extrapolated from the lowest two LES grid levels. The resulting nondimensional shear profile is virtually indistinguishable from the orange curve in Fig. 3 obtained using the same TBLE but with a constant HPPGF within the wall FIG. 3. The nondimensional shear in the LES runs of an idealized ABL applying four different approaches to the modeling of near-surface turbulence. The results are averaged over the horizontal direction and over t 5 5-6 h of the simulations. The gray dotted line is the filtered LOTW in (17). layer (not shown). The insensitivity of TBLE to the methods of estimating HPPGF is expected for an inertial-driven ABL where the perturbation pressure gradient is not a leading term in any of the momentum budget equations.

Although the TBLEs slightly compromise rather than improve the reproduction of LOTW, they efficiently remove the physically unrealistic assumption of instant, local equilibrium between unresolved surface shear stress and resolved nearsurface strain rate. Theoretically, an instant, local equilibrium between the near-surface wind and the surface shear stress implies that their relationship is characterized by a scalar (e.g., a drag coefficient), inconsistent with the physical expectation of a tensor (Wyngaard 2004). Figures 4b and4c shows that the instantaneous, local stress-strain-rate relationship computed using the TBLEs (red cross markers) can deviate remarkably from those obtained using the semi-slip boundary condition (blue dots). Similar deviations from the equilibrium relationship are reported by the filtered direct-numericalsimulation (DNS) results of a smooth-wall boundary layer (see Fig. 8b in Yang et al. 2019). Note that such remarkable deviation is not captured by the turbulence-memory model (red cross markers in Fig. 4a).

This test confirms that the TBLEs are implemented into CM1 LES correctly. In spite of no improvement in reproducing the LOTW, Fig. 4 demonstrates a major advantage of using the TBLE approaches while keeping a consistent result for the LOTW compared to the semi-slip and turbulence memory approaches (Fig. 3).

Applying TBLEs coupled with CM1 LES to simulating an idealized tornado

Theoretically, a semi-slip boundary condition assuming horizontal homogeneity is not realistic for a highly inhomogeneous flow like a tornado. Testing TBLEs on the LES of an idealized tornado reveals the potential error coming from the semi-slip boundary condition applied to a highly curved flow. Specifically, this section seeks to determine how the turbulencememory model and TBLEs change (i) the relationship between the first-level wind and the surface shear stress and (ii) the intensity of a simulated tornado, compared to that with a semislip boundary condition.

a. Simulation configuration

To avoid the cost of simulating an entire tornadic supercell, the so-called "Fiedler Chamber" approach (Fiedler 1995) is used to obtain idealized tornado-like vortices (hereafter called "idealized tornadoes") in a controlled manner that enables comparison of the various near-surface turbulence models. In this approach, an updraft is driven by a constant upward body force in a rotating domain, which results in the formation of a tornado (see Fig. 2 of Rotunno et al. 2016). Such simulation is characterized by a swirl ratio:

where V is the domain's rotating rate, l r is the horizontal scale of the updraft forcing region, and W is the updraft velocity scale (see more details in Rotunno et al. 2016). Here a value of S r 5 0.01 is taken. Each simulation is initialized as a neutral atmosphere at rest. To enable development of turbulent motions, random potential temperature perturbations within 60.25 K are applied to the initial and lateral boundary conditions as used by Wang et al. (2020). Each CM1 LES run employs a domain of 24 km 3 24 km 3 15 km with an inner fine mesh of 4 km 3 4 km 3 1 km (located at the bottom center of the outer domain). The fine mesh grid spacing is 10 m in each direction. From the inner-mesh's edges, the horizontal grid spacing is stretched gradually to reach 190 m at the lateral boundaries, while the vertical grid spacing is stretched gradually to reach 190 m at the top boundary. Four models of near-surface turbulence are tested. Similar to section 3a, the TBLEs embed 20 levels in the wall layers and derive surface shear stress using u and y at the fifth TBLE grid level (i.e., z 5 z wl ≃ 1.25 m), with a roughness length set to z 0 5 0.2 m (the same as the value used by Wang et al. 2020). The simulations using the semislip boundary condition and the turbulence-memory model are the same as cases "OLD-DX10" and "NEW-DX10-g0.1" in Wang et al. (2020), respectively, except that the momentum advection no longer employs a WENO scheme here (Wang et al. 2021). Each simulation runs for 6 h, with an adaptive time step computed to maintain numerical stability. All simulations reach a quasi-steady state after t 5 4 h.

b. Examining the TBLE assumptions about the HPPGF and subgrid-scale mixing

Analyses similar to those in section 3b are performed to examine the TBLE assumptions on the HPPGF. To focus on tornado dynamics, statistics are computed for only a 600 m 3 600 m area in the center of the horizontal domain. Figures 5a-f shows that the leading terms in the w budget are the VPPGF, WTEND, and HADV, meaning that (3) is invalid. Because estimating the variation of WTEND is nontrivial, the p sfc is estimated using the extrapolation (13) rather than integrating the w budget vertically. Although (3) is invalid, Fig. 5g shows that 83.1% of DHPPGF/HPPGF 1 values are below 0.05% and 91.1% of DHPPGF/HPPGF values are below 0.1, suggesting that the assumption of a constant HPPGF with height within the wall layer is approximately valid in most locations. Such results are consistent with the interpretation that (3) is a sufficient FIG. 5. As in Fig. 2, but from the simulation of an idealized tornado within a 600 m 3 600 m area at the domain center. Note that the plotted distributions do not include a few grid points with VPPGF values much larger than 10 m s 22 , which are balanced by similarly large negative values for HADV and/or VADV.

but not necessary condition for assuming a constant HPPGF with height.

To examine how the w budget varies with the distance from the tornado center, Figs. 6a and6b shows the timeazimuthal-mean w budget with respect to the tornado center identified using the instantaneous maximum vertical vorticity at z 5 0.5Dz. The VPPGF is mainly balanced by vertical advection in the region very close to the tornado center (r # 50 m, where r is radius), while HADV becomes increasingly important as r increases (see the purple solid line and blue dashed line in Fig. 6b). At r , 100 m, the timeazimuthal-mean values of DHPPGF/HPPGF 1 are above 0.1 (Fig. 6c), where VPPGF is strong, and the vertical velocity is comparable to the horizontal velocity down to z 5 1.5Dz (Fig. 6d). These results suggest that assuming constant HPPGF with height within the wall layer becomes invalid near the tornado center.

In addition to assuming constant HPPGF within the wall layer, the assumption of neglected horizontal mixing suggested by ( 2) also needs to be examined near the tornado center. Figure 7 shows the time-azimuthal-mean horizontal velocity budget where each term is calculated as

where Term u and Term y are the corresponding components of the u and y budgets, respectively. The horizontal subgridscale mixing (red line in Fig. 7) can be ignored compared to the vertical one (green line in Fig. 7) as assumed by ( 2) roughly when r . 150 m. In the region where r # 150 m, ignoring horizontal subgrid-scale mixing may lead to underpredicted stress. In summary, the TBLE assumptions about the HPPGF and subgrid-scale mixing appear to be valid for the idealized FIG. 6. (a) time-azimuthal-mean (with respect to the tornado's center) w budget at z 5 Dz and (b) each term scaled by the VPPGF term. The abbreviation of each term's name is as in Fig. 2. (c) The difference between the surface and first-grid-level HPPGH (blue line, see the definition in the caption of Fig. 2g, and the gray dashed line denotes the value of 0.1) along with the first-level HPPGF (red line). (d) The velocity's aspect ratio at the lowest three levels within LES domain. The results are all averaged over t 5 4-6 h. tornado simulation except for the region close to the tornado center where the flow is not mainly horizontal and the horizontal shear is comparable to the vertical shear.

c. Results

Figure 8 shows that simulations employing the four models of the near-surface turbulence have similar temporal-movingaveraged maximum updrafts and horizontal velocities once the idealized tornado reaches a quasi-steady state. Before reaching the quasi-steady state, the tornado produced by the Cabot TBLE simulation is apparently stronger than the other three simulations (comparing the orange lines to the other lines in Fig. 8). The sensitivity of TBLE results to the eddy-viscosity model reveals that a tornado's intensity can be significantly changed simply by using a different eddy-viscosity model (even with the same near-surface-turbulence model). Without the temporal moving average (Fig. 9), both TBLE simulations display higher instantaneous updrafts and horizontal velocities (see the spikes of orange and blue lines in Fig. 9).

The instantaneous vertical velocity fields (Fig. 10) show that the nonequilibrium models of near-surface turbulence (including the TBLEs and the turbulence-memory model) yield more turbulent near-surface wind fields than the equilibrium model (i.e., the semi-slip boundary condition). As explained by Wang et al. (2020), the turbulence-memory model captures an additional pathway of dynamic instability development, which enhances energy cascade and consequently the intensity of resolved turbulence. In the runs employing FIG. 7. As in Fig. 6a, but for the horizontal velocity budget at z 5 Dz/2 and with a log-scale y axis and absolute values, where UTEND is the tendency of horizontal velocity calculated as the residual of summing all other terms. TBLEs, the enhanced resolved turbulence near the surface can be explained by the backscatter of TKE (discussed later).

To reveal the quasi-steady-state stress-strain-rate relationship for different models of near-surface turbulence, Fig. 11 shows time-averaged horizontal wind at the lowest LES grid level and the surface shear stress with respect to the tornado center. Note that the tornado center is not necessarily the domain center because the tornado drifts slightly. When the directions of horizontal wind and surface shear stress are described, "radial" and "tangential" directions of a cylindrical coordinate system are used with respect to the tornado center, while "outward" implies having a positive radial component (i.e., away from the tornado center) and "inward" implies having a negative radial component (i.e., toward the tornado center). All simulations yield similar wind directions (represented by blue arrows in Fig. 11): the velocity points inward at r 5 50 m but is almost purely tangential at r 5 20 m. All simulations yield similar surface-shear-stress directions at r 5 50 m, but substantially different surface-shear-stress directions arise at r 5 20 m (represented by yellow arrows in Fig. 11). Specifically, the surface shear stress given by a semislip boundary condition is directed opposite the horizontal wind at the lowest LES grid level (Fig. 11a) as constrained by the instant, local equilibrium. The surface shear stress given by the turbulence-memory model has a component contributing to an air parcel's centripetal acceleration (Fig. 11b) owing to the nonzero lifetime of turbulent motions (Wang et al. 2020). Both TBLE simulations produce surface shear stress directed outward (Figs. 11c,d). This is because, while the horizontal wind is nearly tangential at the lowest LES grid level (Figs. 12a,b), it is directed inward within the wall layer (Figs. 12c,d) and the shear stress is directed opposite the wall layer velocity.

To explore the cause of the short-lived but intense updrafts (i.e., the spikes in Fig. 9) produced by CM1 LES employing TBLEs, the temporal average of each instance when the maximum updraft exceeds 180 m s 21 is examined. Figure 13 reveals that at the moment of an intense updraft, the horizontal wind at the lowest LES grid level points inward at r 5 20 m (cf. Figs. 11c,d) and the surface shear stress has a component contributing to an air parcel's centripetal acceleration, similar to that given by turbulence memory. The inward-directed velocity at r 5 20 m is associated with a maximum updraft at the center at the lowest LES level (Figs. 14a,b,cf. Figs. 12a,b), indicating the existence of an end-wall vortex (Church et al. 1977;Fiedler and Rotunno 1986;Lewellen and Lewellen 2007;Rotunno 2013;Davies-Jones 2015). The fifth TBLE grid level shows a ring of downdraft (represented by the dark blue regions in Figs. 14c,d) with the horizontal winds less inward-directed than those at the lowest LES grid level (comparing the white arrows in Figs. 14c,d to those in Figs. 14a,b), thus yielding a surface shear stress contributing to an air parcel's centripetal acceleration near the tornado center.

Because Cabot TBLE behaves most differently from the semi-slip boundary condition, the simulation using Cabot TBLE is hereafter used to explore the mechanism of the briefly intensified tornado. To further explore the moment at which the surface shear stress contributes to an air parcel's centripetal acceleration, time-azimuthal-mean flows within the LES domain and the wall layer are presented in Figs. 15a and15b. The LES flow field shows an end-wall vortex (Fig. 15a), while the wall-layer flow shows a corner backflow at r , 50 m and z , 2 m (Fig. 15b). The corner backflow is produced by a near-surface outward-directed HPPGF (Fig. 15c; see the blue line with negative values at 30 , r , 50 m) which apparently dominates the inward advection term in this region of the wall layer. Contrary to the average of intense-updraft moments, the quasi-steady-state tornado never shows outward directed HPPGF near the surface (see the blue line in Fig. 16c), so FIG. 9. As in Fig. 8, but the values are the original output every second instead of the 5-min moving average.

the corner backflow and end-wall vortex are not observed in the quasi-steady state (Figs. 16a,b). Further inspection of the outward directed HPPGF, which is responsible for developing the wall-layer corner flow identifies a local perturbation pressure maximum at r ≃ 30 m and a local perturbation pressure minimum at r ≃ 50 m (see the gray line in Fig. 15c). The local perturbation pressure maximum at r ≃ 30 m coincides with strong convergence. Specifically, the inward directed HPPGF outweighs the centrifugal force for r . 60 m, accelerating air parcels toward the tornado center. At r , 60 m, the centrifugal force outweighs the inward directed HPPGF, decelerating the inward moving air parcels. As a result, strong horizontal convergence occurs at r , 60 m, which induces large pressure perturbations (a phenomenon reported by Ward 1972, owing to a physical mechanism explained by Rotunno and Klemp 1982). A positive feedback mechanism between the pressure perturbations and the horizontal convergence leads to rapid intensification of the tornado (note that the sum of the budget terms is not small, consistent with an unsteady state). The quasi-steady state also shows a slight increase in perturbation pressure (Fig. 16c, see the gray line's flattened slope within 30 , r , 60 m), but the increase is not as large as that associated with the end-wall vortex.

In addition to the corner backflow, the dislocation of the surface shear stress' spiral center (defined by the maximum vertical vorticity of the wall-layer flow where the surface shear stress is derived) with respect to the tornado center (Fig. 17) also contributes to the intense updrafts. Specifically, the air parcel between the two centers experiences surface shear stress that contributes to the centripetal acceleration (Fig. 17b) and approaches the tornado center more easily. About 10 s after the dislocation of the two centers, an end-wall vortex with a strong near-surface updraft develops (Fig. 17d), and the centers begin to realign (Fig. 17e). About 5 s after the occurrence of the end-wall vortex, the two centers become collocated (Fig. 17h), and the end-wall vortex is replaced with a strong near-surface downdraft (Fig. 17g).

When the angles between horizontal winds at the lowest LES grid level and the surface shear stresses are smaller than 908 (between the two centers in Figs. 17b,e,h), a backscatter of kinetic energy from unresolved to resolved scales takes place (see the observation and explanation in Sullivan et al. 2003). The strength of backscatter can be quantified approximately using the dot product between surface shear stress and shear at the lowest LES grid level (i.e., 2T sfc derived in the appendix). Figures 17c,17f, and 17i shows that the backscatter peaks when the two centers are most apart from each other and decreases as the two centers become close. The semi-slip boundary condition always yields a positive-definite T sfc , meaning that the near-surface turbulence always removes kinetic energy from the system. The turbulence-memory model introduces an additional surface-shear-stress component normal to the horizontal wind at the lowest LES grid level (Wang et al. 2020), which does no work on an air parcel and therefore does not influence T sfc . The TBLEs allow even more flexible surface-shear-stress directions, and consequently no longer require T sfc to be positive-definite, thus permitting the backscatter.

The abovementioned two mechanisms (the wall-layer corner backflow and the dislocation of the two centers) are further investigated by their temporal correlation with the maximum updraft. The intensity of the wall-layer corner flows is measured using the maximum outward directed radial flow in the wall layer, and the dislocation is measured using the distance between the two centers. Figure 18 reveals that the wall-layer corner backflow's intensity is highly correlated with the maximum updraft (Fig. 18a), and the high correlation occurs mainly during the transient stage (blue circle in Fig. 18a, cf. orange line in Fig. 11b). On the other hand, the dislocation of the two centers is uncorrelated with the maximum updraft (Fig. 18b), although the causality cannot be fully excluded. The correlation may be affected by that with the corner backflow (i.e., the dominant mechanism), and there is a time lag from the misalignment of two centers to the development of the end-wall vortex (Fig. 17), thus reducing the representativeness of the instantaneous correlationfoot_7 ; a possible reason is the time lag of the development of an end-wall vortex to the two centers' dislocation (;10 s as revealed by Fig. 17), which cannot be represented by the instantaneous correlation.

Conclusions

This work investigates the potential consequences of employing an assumption of statistically steady and horizontally homogeneous near-surface turbulence (i.e., the equilibrium approach) in LES of unsteady and horizontally heterogeneous atmospheric flows. CM1 LES runs are performed using four different models of the near-surface turbulence: (i) the semi-slip boundary condition based on MOST (which is an equilibrium approach), (ii) a model recently proposed to account for the influence of turbulence memory associated with curved trajectories (which is no longer an equilibrium approach), (iii) a nonequilibrium approach, TBLE, adopted from the engineering community (Balaras et al. 1996), and (iv) another TBLE, but with a different eddy-viscosity model for closing the RANS equations (Cabot 1996). As an approach originally proposed for LES of smooth-wall boundary layers, the TBLE is modified to become applicable to atmospheric flows above rough surfaces when it is implemented into CM1.

The sensitivity tests of the four models of near-surface turbulence are conducted first for an idealized ABL, to validate the TBLE implementation, and then for an idealized tornado characterized by rapid acceleration, strongly curved air parcel trajectories, and substantial horizontal heterogeneities. For an idealized ABL, the turbulence-memory model behaves approximately the same as an equilibrium model owing to little curvature of air parcels' trajectories. The TBLEs yield mean-wind profiles similar to those obtained using the other two models, with a slightly worse reproduction of the theoretical mean-wind profile right above the wall layer. Nevertheless, TBLEs effectively break the instant, local equilibrium between the surface shear stress and the horizontal wind at the lowest LES grid level, yielding results qualitatively consistent with filtered DNS results (e.g., Yang et al. 2019). For an idealized tornado, TBLEs yield much stronger instantaneous velocity maximum values than the semi-slip boundary condition and the turbulence-memory model. An average across all intense-updraft events (maximum w . 180 m s 21 ) shows that the surface shear stress computed using TBLEs contributes to an air parcel's centripetal acceleration, like what occurs with the turbulence-memory model. The local outward-directed HPPGF creates a corner backflow (i.e., positive radial velocity) within the wall layer, whose intensity is correlated with the maximum updraft. The dislocation of the surface shear stress' spiral center from the tornado center permits backscatter of kinetic energy from unresolved to resolved scales, thus increasing the resolved turbulence near the surface.

Whether the TBLE outperforms the semi-slip boundary condition in reproducing tornado dynamics remains unknown, because current observational technology cannot measure near-surface turbulence in a tornado. The lack of observational data has been a main reason why LES remains a critical tool in understanding convective storms. The current results show that by replacing an equilibrium approach with solving the fluid dynamics equations (i.e., the TBLE), one can obtain a substantially different tornado. Furthermore, even a small change in the eddy-viscosity model employed by the TBLE can lead to large variability of the simulated tornado. Note that the TBLE has not completely eliminated all assumptions implicitly taken by an equilibrium approach. For example, the assumption of a constant HPPGF with height within the wall layer is not uniquely taken by TBLE, but also taken by equilibrium approaches through the assumption of statistically steady and horizontally homogeneous nearsurface turbulence (see details of derivations in Pope 2000, chapter 7.2.1). TBLEs have been used by a number of studies without challenging the assumption of a constant HPPGF with height (e.g., Balaras et al. 1996;Cabot 1996;Cabot and Moin 2000;Wang and Moin 2002). Park and Moin (2014) evaluated the vertical pressure gradient but found it small in their case (personal communication; not shown in their paper).

To our knowledge, this work is the first study that emphasizes the necessity of examining the assumption of a constant HPPGF with height for modeling near-surface turbulence. Other example assumptions include an approximately horizontal velocity at the lowest LES grid level, the negligible horizontal subgrid-scale mixing, a characteristic length scale depending only on the distance from the surface, and a single parameter z 0 to represent the influence of canopy roughness. In other words, the spread of simulated tornado results obtained in this work may represent only part of the deviation between LES results computed using the semi-slip condition and reality, especially near the center of the tornado where some assumptions are obviously invalid. Nevertheless, because these assumptions have become explicit in the TBLE, different approaches of estimating the associated quantities can be tested (e.g., different methods of estimating HPPGF within the wall layer have been tested in section 3c). As a first attempt to apply TBLE to atmospheric flows, this work has not considered the influence of temperature stratification on near-surface turbulence. Restoring the internal energy equation to the wall-layer RANS equations needs to be the first step toward extending the TBLE application to true tornadoes within parent storms. The closure schemes for the RANS equations, including the eddy-viscosity model and the calculation of surface shear stress, also need to account for the influence of temperature stratification.