## Abstract

The impact of atmospheric boundary layer (ABL) interactions with large-scale stably stratified flow over an isolated, two-dimensional hill is investigated using turbulence-resolving large-eddy simulations. The onset of internal gravity wave breaking and leeside flow response regimes of trapped lee waves and nonlinear breakdown (or hydraulic-jump-like state) as they depend on the classical inverse Froude number, Fr^{−1} = *Nh*/*U*_{g}, is explored in detail. Here, *N* is the Brunt–Väisälä frequency, *h* is the hill height, and *U*_{g} is the geostrophic wind. The results here demonstrate that the presence of a turbulent ABL influences mountain wave (MW) development in critical aspects, such as dissipation of trapped lee waves and amplified stagnation zone turbulence through Kelvin–Helmholtz instability. It is shown that the nature of interactions between the large-scale flow and the ABL is better characterized by a proposed inverse compensated Froude number, = *N*(*h* − *z*_{i})/*U*_{g}, where *z*_{i} is the ABL height. In addition, it is found that the onset of the nonlinear-breakdown regime, ≈ 1.0, is initiated when the vertical wavelength becomes comparable to the sufficiently energetic scales of turbulence in the stagnation zone and ABL, yielding an abrupt change in leeside flow response. Finally, energy spectra are presented in the context of MW flows, supporting the existence of a clear transition in leeside flow response, and illustrating two distinct energy distribution states for the trapped-lee-wave and the nonlinear-breakdown regimes.

## 1. Introduction

A significant body of work exists in the literature on mesoscale stably stratified flows over topographic features, including analytical, experimental, and numerical efforts (e.g., Queney 1948; Lilly 1978; Klemp and Lilly 1978; Durran 1986; Durran and Klemp 1987; Smolarkiewicz and Rotunno 1989, 1990; Baines and Hoinka 1985; Baines 1998; Kim et al. 2003). However, the vast majority of studies on this topic have focused on atmospheric conditions at elevations significantly above hill height. The motivation was often twofold, with the emphasis being either on aviation safety, as mountain waves (MW) are one possible catalyst in aircraft accidents (e.g., Lilly 1978; Sharman et al. 2012), or on the broader effects of MW drag on numerical weather prediction and global-scale atmospheric dynamics (e.g., Broad 1996; Welch et al. 2001; Xu et al. 2013). In addition, the added complexity of atmospheric boundary layer (ABL) turbulence from both a theoretical and numerical standpoint likely precluded explicit inclusion of ABL interactions. As such, the influence of the ABL on MW development and dynamics was neglected by omission or via overly simplified boundary layer parameterization schemes. Nonetheless, these and other similar significant efforts have resulted in the foundational body of work in characterization of MW regimes via the inverse Froude number, or nondimensional hill height, Fr^{−1} = *Nh*/*U*_{g}, where *N* is the Brunt–Väisälä frequency, *h* is the hill height, and *U*_{g} is the geostrophic wind.

It has been suggested in numerical studies by Richard et al. (1989), Ólafsson and Bougeault (1997), and Peng and Thompson (2003) that the influence of the ABL is to reduce MW development and inhibit wave breaking aloft. In Ólafsson and Bougeault (1997), surface friction was shown to suppress MW, and rotation would either suppress MW (Fr^{−1} < 1.4) or enhance MW for larger nondimensional hill heights. A boundary layer parameterization was used in Peng and Thompson (2003) to show that ABL height substantially impacts MW development. Interaction between trapped lee waves and the ABL was further studied with a mesoscale model by Jiang et al. (2006). A two-dimensional theory for idealized boundary layer effects was developed in Smith et al. (2006), followed by a three-dimensional linear theory in Smith (2007). Each of these studies, without explicitly accounting for any turbulent processes, contributes to the general consensus that the ABL inhibits MW activity. Epifanio and Qian (2008) in a numerical study explored the dynamics of turbulence development due to wave breaking and concluded that the mean flow is rapidly dissipated primarily because of generation of this turbulence. However, these results are independently inconclusive, given (i) the coarse resolutions used (250–2500 m) and (ii) the lack of resolved initial and inflow turbulence results in a misrepresentative flow solution (see section 4). In addition, Doyle et al. (2011) show that a free-slip surface boundary condition as used in Epifanio and Qian (2008) in the absence of a surface drag parameterization scheme fails to produce any resolved small-scale turbulence, resulting in nearly identical results to the typical mesoscale modeling solution.

With increasing computational capabilities, large-eddy simulation (LES) has provided a higher-fidelity tool for studying ABL influence on MW development. Smith and Skyllingstad (2009) used LES to study upstream boundary layer stability influence on MW development. The authors showed that a convective (unstable) boundary layer inhibits the onset of MWs, while an upstream stable boundary layer leads to stronger downslope winds and rotor intensification, providing insight into the role of surface heat flux on MW regime development. Then Smith and Skyllingstad (2011) showed that the effect of a strong noncapping inversion exacerbates MW development with decreasing inversion height. These two efforts not only lend credibility to the applicability and relevance of LES modeling as a tool in MW studies but also make a pioneering contribution through demonstration of turbulent-inflow ABL influence on MW development.

This manuscript aims to elucidate the fundamental fluid dynamics of ABL interactions with stratified flow over an idealized hill. Specifically, we examine both the onset of internal gravity wave (IGW) breaking and leeside flow response as a function of Froude number. Additionally, we explore in detail how the inclusion of a turbulent ABL supplied as both initial and boundary conditions affects classical MW regimes. Finally, we discuss new facets of the role of turbulence in leeside flow response through the analysis of energy spectra, a classical turbulence quantification metric that has not been explored in the existing MW literature. The remainder of this manuscript is organized as follows. Section 2 provides details of model formulation and simulation setup. Section 3 presents detailed analysis of LES-resolved velocity and turbulent kinetic energy (TKE) fields in addition to vertical profiles of TKE and key TKE budget terms. Section 4 proposes a compensated Froude number, which better predicts both the onset of IGW breaking and leeside flow response and demonstrates the importance of ABL-inclusive LES simulations. TKE spectra and flow structures are presented in section 5. Finally, section 6 is devoted to the conclusions of this work and important considerations for future efforts.

## 2. Model formulation and simulation setup

Simulations herein were performed using the modernized version of the High Gradient (HiGrad) large-eddy simulation model (Sauer 2013). The subsections here describe first the LES model formulation in HiGrad and then the specific simulation setup used in this study to produce MW flows.

### a. Model formulation

Fundamentally, LES models solve the filtered Navier–Stokes equations in a manner such that flow features (or eddies) containing most of the energy of the flow are explicitly resolved by the discrete grid. Critically, all LES modeling efforts must provide sufficient resolution to ensure resolved flow into the inertial subrange where eddies are expected to follow the −5/3 slope of Kolmogorov theory for the energy cascade of three-dimensional isotropic turbulence (Kolmogorov 1941). The influence of less-energetic, subgrid-scale flow features is captured in HiGrad via a gradient approximation model for the subgrid-scale stress tensor. Through inclusion of a prognostic equation for subgrid turbulent kinetic energy, a variable eddy viscosity is applied to capture the mixing effects from flow fluctuations below the chosen grid resolution (Lilly 1966, 1967).

The conservation equations in HiGrad are implicitly filtered at the grid-length scale by the numerical discretization, resulting in the following flux conservative form of the filtered compressible LES equations for mass, momentum, and energy (potential temperature), where denotes the grid-filtered value of the prognostic variable *x*:

Here, is air density, and is the component of velocity in the *i*, *j* = 1, 2, 3 directions, indicating zonal (streamwise), meridional (spanwise), and vertical, respectively. A generalized source term for momentum is denoted by . The energy equation [Eq. (3)] is posed in terms of potential temperature , defined by , where and are temperature and pressure of the control volume. The surface pressure is denoted by , *R* is the ideal gas constant, and is the specific heat of dry air at constant pressure. In Eq. (2), the perturbation pressure and density are denoted by and , respectively, and indicate the departure from the hydrostatically balanced environmental base state. The subgrid-scale Reynolds stress tensor is parameterized following Lilly (1966, 1967) for momentum and Deardorff (1980) for heat:

where is the dynamic eddy viscosity, *K* is the total subgrid-scale TKE, and Pr is the turbulent Prandtl number. The model constant is set to , following Moeng and Wyngaard (1988). The subgrid scale or filter length includes a modification from the nominal grid size, , in order to take into account possible reductions of as a result of thermal stratification (Deardorff 1980),

which is in turn used to account for stability effects in the turbulent Prandtl number definition, . The prognostic equation for *K* is given by

where the left-hand side represents the time rate of change and advective tendency of unresolved TKE. The terms on the right-hand side of this equation are, from left to right, shear generation, buoyancy source/sink, turbulent transport, and dissipation. Turbulent transport is modeled via downgradient diffusion, and the dissipation term is based on the Kolmogorov hypothesis (for more details, see Deardorff 1980; Moeng 1984; Muñoz-Esparza et al. 2014a).

### b. Simulation setup

We employ a one-way nested approach to perform simulations of flow over an idealized hill with a fully developed turbulent ABL in both initial and evolving boundary conditions. A flat-terrain outer domain (parent) is used to produce developed turbulence through periodic lateral boundary conditions and concurrently imposed as streamwise and top boundary conditions for each prognostic variable of the nested domain, which is periodic in the spanwise direction to accommodate the two-dimensional hill setup (see Fig. 1). For all the simulations performed herein, we use a uniform-in-height and constant-in-time geostrophic wind aligned in the *x* direction of (*U*_{g}, *V*_{g}) = (0, 0) m s^{−1}.

The surface boundary of the nested domain contains a two-dimensional hill defined by

where *h* is the mountain height (see Table 1), is the location of the mountain center, and the half-width of *a* = 3000 m. Surface turbulent fluxes are parameterized using the Monin–Obukhov similarity theory (Monin and Obukhov 1954), which relates momentum fluxes at the surface to velocity components at the first grid point above the ground in the following manner:

where *i* = 1, 2 denotes the zonal and meridional (*x* and *y* direction, respectively) components of velocity. The second index, 1, designates the first vertical grid point; is the momentum exchange coefficient given by for neutral conditions; and *κ* is the von Kármán constant (=0.4). This formulation does not include the stability correction term since we focus on a neutrally stratified ABL, where the surface heat flux is zero. We set the surface roughness length to 0.1 m in all the cases, representative of the unimproved landscape typically found over mountainous terrain (see Stull 1988).

All simulations have a neutrally stratified boundary layer near the ground of 300 K over an absolute depth of m (i.e., neglecting the hill). Initial conditions for the velocity field components, however, are interpolated from the parent (flat terrain) domain to the nested (hill inclusive) domain with reference to above-ground-level height. Taken in combination, this pair of initial conditions allowed the turbulent (shear only) boundary layer to be initially draped over the entire domain, including the hill, but maintained the initial horizontal homogeneity of potential temperature with height such that the presence of the hill could disrupt the ambient stratification and instigate MW activity. Different Froude regimes were prescribed by modifying the strength of a constant-with-height potential temperature gradient for *z* > *z*_{i} (see Fig. 2), and hence the Froude number as indicated in Table 1. In this manner, we explore a large range of the parameter space Fr^{−1} = 0.9–3.4, albeit through a set of eight discrete values. Table 1 includes values of inverse compensated Froude number, , which is introduced in section 3. The dimensions of the parent domain are × × = 36 km × 6 km × 16 km, while the nested domain is slightly reduced: × × = 32 km × 4.8 km × 12 km. The horizontal grid spacing is 20 m in both the *x* and *y* directions. The grid is stretched vertically from 20 m at the surface to 530 m at the domain top following a cubic polynomial relationship.

Typically, LESs of ABLs using periodic lateral boundary conditions are performed imposing a large-scale pressure gradient force on the momentum equation. This force can be parameterized as a perturbation from geostrophic conditions that represents the veering of the horizontal velocities with height as a result of its balance with the Coriolis force (e.g., Moeng 1984; Kosović and Curry 2000; Muñoz-Esparza et al. 2014b). The idealized two-dimensional flow configuration we are interested in reproducing herein requires the flow to be perpendicular to the two-dimensional hill, hence precluding the application of such Coriolis parameterization of the large-scale pressure gradient force. Instead, we use a similar method to that proposed by Ólafsson and Bougeault (1997) to impose a large-scale pressure gradient force that results in a realistic wind speed profile in the streamwise direction at any height within the boundary layer and remains constant with height above. This method is based on employing the projection of balanced large-scale pressure gradient and Coriolis forces onto the streamwise momentum profile at each discretized height in the domain. This forcing term is given by

In Eq. (10), ∂*p*/∂*x*, ∂*p*/∂*y*, and ∂*p*/∂*z* are the corresponding large-scale pressure gradients in the *x*, *y*, and *z* directions, and and are the equilibrium horizontal velocities obtained from a reference periodic turbulent ABL simulation using Coriolis parameter *f* and a geostrophic wind *U*_{g} aligned in the *x* direction. Here, *f* = 1.181 × 10^{−4} s^{−1} corresponding to a latitude of 54°N is used. By applying the pressure gradient term described above, a logarithmic profile boundary layer aligned in the *x* direction is obtained within the parent domain such that resolved TKE is decreasing with height once the influence of the parameterized surface drag disappears (*z* > 50 m; see Fig. 2). Figure 2 provides instantaneous profiles (gray) and mean profiles (blue) at an arbitrary location of the parent domain for streamwise velocity, resolved TKE, and potential temperature representative of initial and boundary conditions used in the nested domain.

The governing equations described in the previous section are discretized using the finite-volume technique under the third-order accurate-in-space quadratic upstream interpolation for convective kinematics (QUICK) advection scheme (Leonard 1980) and third-order accurate Runge–Kutta time integration scheme (Durran 2013). A Rayleigh-damping layer was applied in the parent domain in the uppermost region to avoid wave reflection from the domain top. As a result of HiGrad being a fully compressible code, the time step was chosen to be *dt* = 0.05 s in order to satisfy *dt* < min(*dx*, *dy*, *dz*)/*c*, where *c* is the speed of sound. The parent domain was run for 10 h of simulated time in order to establish a fully developed turbulent boundary layer, and then the nested domain was initiated. Both domains were run concurrently for an additional ≈2.5 h. This time was found to be sufficient for the MW flows to reach a quasi-steady state. Instantaneous results shown in the rest of the manuscript correspond to *t* = 2.5 h or nondimensional time . All time-averaged fields and profiles are taken over the last 10-min interval: .

## 3. LES-resolved flow dynamics

The benefit of utilizing LES for studying interactions between the ABL and large-scale stratified flow over a terrain feature lies in the ability of LES to explicitly resolve the majority of turbulent flow dynamics. By eliminating parameterizations (or gross approximations) of complex turbulent processes, the dynamic interactions between the ABL and larger-scale flow over a terrain feature leading to the development of enhanced shooting flow, IGW breaking, and leeside flow response are not constrained by assumptions made in the parameterization of turbulence. Through LES, a less approximated and therefore more informative study can be accomplished in order to better understand the influence of a turbulent ABL on these intricate features of MW flows.

### a. Streamwise velocity fields

In Fig. 3, instantaneous streamwise velocity fields are plotted as filled contours in the spanwise–centerline *x*–*z* plane for six different Froude numbers in the range Fr^{−1} = 0.9–2.7 and hill height *h* = 900 m (see Table 1). For the weakest stratification case corresponding to Fr^{−1} = 0.9, there is only a mild enhancement of zonal velocity above the hill and absolutely no development of stagnation zone flow reversal or subsequent IGW breaking. On the lee side of the hill, the flow is nearly undisturbed with no leeside rotor development. For the Fr^{−1} = 1.3 case, an enhanced shooting flow roughly twice the geostrophic wind speed develops just beyond the hill crest, and stagnation zone flow reversal generates turbulence because of IGW breaking. Multiple rotors (regions of flow reversal in the turbulent ABL) develop in the leeside flow. The shooting flow interaction with the stagnation zone above and the ABL below produces a trapped lee wave that is weakly dissipated with increasing downstream distance.

The Fr^{−1} = 1.6 and 1.8 cases in Fig. 3 exhibit stronger enhancement of shooting flow to approximately 3 times the geostrophic wind speed and clear regions of stagnation zone reversed flow, leading to overturning of potential temperature contours and IGW breaking, with second-level IGW breaking in the latter case. The flow field on the lee side of the hill in both cases develops multiple rotors, with the first leeside rotor containing the highest reversed winds. These results show that the resultant trapped lee wave is the downstream manifestation of the shooting flow trapped between the turbulent ABL and the critical layer of Kelvin–Helmholtz (KH)-induced and downstream advected turbulence of the stagnation zone.

The stronger stratification cases, Fr^{−1} = 2.3 and 2.7 in Fig. 3, show the maximum velocity of the shooting flow is greater than 3 times the geostrophic wind speed, and multiple levels of stagnation zone reversed flow result in IGW breaking. In both cases, only a single rotor structure is present, and the stagnation zone and ABL layers merge into a single, deep layer of turbulent flow on the downstream side of the single rotor. For clarification, throughout the remainder of this paper we denote a leeside nonlinear-breakdown response scenario as one in which only a single rotor is present, whereas the presence of multiple rotors is our criteria for cases we term trapped-lee-wave scenarios.

With increasing Fr^{−1}, either the horizontal length scale of trapped lee waves or the single leeside rotor decreases, as the upstream edge (or separation point) of the first rotor is pushed downstream. The growing enhancement of the shooting flow more readily pushes the leading edge of the turbulent ABL farther down the lee side of the hill. As the turbulent ABL resists the shooting flow, this high-speed, essentially laminar wind is lifted off of the surface at the separation point as a result of the local pressure gradient. In addition, the lifted region of the shooting flow carries significant energy toward the critical layer above. The shooting flow can then be reflected back toward the ground by the critical layer, establishing the downstream, descending portion of the rotor.

To verify that the LES approach accurately predicts large-scale MW flow features in addition to capturing the largest and most energetic production and inertial range scales of turbulence in a resolved manner, we consider the vertical wavelength of internal gravity waves *λ*, defined as

where is the Scorer parameter (Scorer 1949), defined as

and is the second derivative of geostrophic wind with respect to the vertical coordinate (in the cases here above the ABL). Visual inspection of the cases shown in Fig. 3 reveals good agreement with the range of predicted vertical wavelengths from λ = 6.2 to 2.1 km over increasing Fr^{−1}, demonstrating that the LES approach, as expected, captures this fundamental large-scale flow feature.

### b. Turbulent kinetic energy fields

Instantaneous spanwise-averaged, resolved turbulent kinetic energy is shown in Fig. 4 as filled contours in the *x–z* plane, providing additional insight into ABL influence on MW development and leeside flow regimes, corresponding to the same cases shown in Fig. 3. The resolved turbulent kinetic energy is given by

where perturbation velocities are calculated as

and is the spanwise (or *y*-direction) mean. The TKE fields shown in Fig. 4 are the spanwise average, , illustrating the coherent structure of TKE that develops in both the stagnation zone and leeside ABL.

In all cases, the mean TKE of the incoming turbulent ABL is on the order of unity, as seen in Figs. 2 and 4. For the weakest stability case where Fr^{−1} = 0.9, no IGW breaking occurs, but the slight increase in streamwise velocity (shooting flow) over the ridge excites leeside TKE levels in the ABL through increased shear generation. In this weakest stability case, the upstream ABL is advected over the hill, with only a slight compression of the upstream vertical profile of TKE over the ridgetop. At *x* ≈ −2.0 km, the TKE occupies a shallower height above ground level, since the descending shooting flow sweeps away ABL TKE and carries in primarily laminar flow from aloft. Farther downstream toward *x* ≈ −1.5 km, the presence of the turbulent ABL imposes a drag on the modest shooting flow, leading to a weak enhancement of downstream TKE in the ABL through shear generation. With increasing stability (increasing Fr^{−1}) the compression of the ABL by the descending shooting flow is more pronounced, resulting in an extremely shallow layer of weak TKE being advected over the ridge. This implies that the role of upstream TKE as advected over the ridge in the onset of IGW breaking or leeside flow response is essentially negligible; however, the mere presence of the upstream turbulent ABL does play an important role in the development of IGW breaking and stagnation zone TKE generation, as discussed more in section 4.

In the top-right panel of Fig. 4, where Fr^{−1} = 1.3, potential temperature contours above the hill indicate overturning leading to buoyancy generation of TKE. Moreover, in this stagnation zone, streamwise velocity is suppressed, while just below, in the shooting-flow region, the streamwise velocity is enhanced. The lower edge of the stagnation zone subsequently shows higher TKE levels through enhanced shear generation since the magnitude of the vertical velocity gradient is increased as the shooting flow region and stagnation zone develop. In this stability case, we again see enhancement of leeside TKE in the ABL, but to a much stronger level (≈20 m^{2} s^{−2}). At the separation point, *x* ≈ −2.5, the shooting flow detaches from the surface and rises back toward the critical layer where substantial TKE exists because of advection from the upstream wave breaking in the stagnation zone and local shear and buoyancy generation. This upper critical layer associated with IGW breaking deflects the shooting flow back toward the ground, contributing to the development of a leeside rotor. Farther downstream at *x* ≈ 5 km, the shooting flow again encounters the surface ABL where roughness enhances drag affects and shear generation of TKE. By this downstream distance, the critical layer and surface layer TKE have combined to absorb the kinetic energy and subsequently decrease the strength of the shooting flow. As the shooting flow again separates from the surface layer at *x* ≈ 5.25 km, a second rotor ensues. The leeside flow regime in this case is clearly that of a trapped lee wave.

For Fr^{−1} = 1.6 and 1.8 in Fig. 4, shooting-flow magnitude increases as vertical wavelength decreases, IGW breaking occurs, and LES-resolved TKE levels show a stronger dissipation of shooting-flow kinetic energy with downstream distance by interaction with the critical layer and surface ABL. Consequently, the trapped-lee-wave regime transitions from weakly dissipative (Fr^{−1} = 1.3) through moderately dissipative (Fr^{−1} = 1.6) to strongly dissipative (Fr^{−1} = 1.8) as stagnation zone and surface ABL turbulence levels increase and more efficiently subsume the increasingly energetic shooting flow.

In the cases of Fr^{−1} = 2.3 and 2.7, the leeside flow response transitions to a state in which a complete nonlinear breakdown of the shooting flow occurs, where only a single rotor exists. In this region of the Froude number parameter space, the stagnation zone turbulence and surface-layer ABL completely subsume the shooting flow. The vertical wavelength of the IGW becomes sufficiently small so that the enhanced stagnation zone and surface-layer ABL turbulence collectively lead to a complete breakdown of the shooting flow. The enhanced magnitude of the shooting flow pushes the separation point farther downstream (*x* ≈ 0 km) and TKE levels in the rotor increase by nearly an order of magnitude. These two stability cases fall into a nonlinear-breakdown state similar, as many other researchers have suggested, to a hydraulic jump (e.g., Long 1953; Durran and Klemp 1987; Doyle and Reynolds 2008).

### c. Profiles of turbulent kinetic energy and budget terms

The vertical profiles of 10-min averaged, resolved TKE at several streamwise locations are presented in Fig. 5 for six different Froude numbers in the range Fr^{−1} = 0.9–2.7 and hill height *h* = 900 m (see Table 1). Columns correspond to four different streamwise locations: *x* = 3, 0, 2, and 4 km from left to right. The top row in Fig. 5 shows the weakest stability case, Fr^{−1} = 0.9, where no IGW breaking occurs. The middle row includes profiles for the weakly, moderately, and strongly dissipative lee-wave flow response scenarios, Fr^{−1} = 1.3, 1.6, and 1.8, respectively, where IGW breaking occurs, but the leeside flow response is such that the shooting flow can persist in the vertical space between the stagnation zone and leeside turbulent ABL. The bottom row in Fig. 5 shows two leeside nonlinear-breakdown cases, Fr^{−1} = 2.3 and 2.7, respectively, where the shooting flow is completely subsumed by the stagnation zone and leeside ABL turbulence. In general, across the six Fr^{−1} values, these downstream locations reasonably correspond to the places where (i) shooting flow dominates near the surface and stagnation zone TKE develops aloft (*x* = −3 km), (ii) boundary layer TKE and stagnation zone TKE begin to interact with the shooting flow (*x* = 0 km), (iii) boundary layer and stagnation zone turbulence have either reached a quasi-steady balance with the shooting flow or entirely merged and subsumed the shooting flow (*x* = 2 km), and (iv) TKE production is steadily diminishing downstream of the leeside slope (*x* = 4 km).

In Fig. 5, the profiles for all six Fr^{−1} at *x* = −3 km (1 km downstream from the peak of the hill) exhibit very little TKE near the surface, where the essentially laminar shooting flow dominates. In the weakest stability scenario, Fr^{−1} = 0.9, there is no development of noticeable levels of TKE in the stagnation zone since there is no IGW breaking. With increasing downstream distance, the intensity and vertical extent of TKE in the leeside ABL initially increases through interaction with the weakly enhanced zonal velocity of the shooting flow, reaching a maximum near *x* = 2 km. By *x* = 4 km, the enhanced zonal velocity has relaxed, and TKE decreases with downstream distance as a result of dissipation and lack of shear generation.

For the weakly dissipative, trapped-lee-wave case of Fr^{−1} = 1.3, the *x* = −3-km profile shows only weak turbulence aloft since this location is at the most upstream edge of the stagnation zone turbulence (see Figs. 3 and 4). For Fr^{−1} = 1.6 and 1.8, the development of TKE in the stagnation zone at *x* = −3 km is apparent with maximum levels of ≈5–10 m^{2} s^{−2}. For these trapped-lee-wave cases, the moderate levels of turbulence in the ABL and stagnation zone along with the vertical wavelength associated with the strength of stratification collectively permit the shooting flow to propagate downstream between the two turbulent regions with only limited transference of KE to TKE through shear mixing. For these three cases, a separation exists between the two highly turbulent regions at each of the downstream locations—*x* = 0, 2, and 4 km—with three inflection points in the vertical profiles of TKE. For the trapped-lee-wave cases, KH-induced TKE in the stagnation zone is lower in magnitude than the enhanced TKE in the turbulent ABL.

In the bottom row of Fig. 5, the leeside nonlinear-breakdown cases of Fr^{−1} = 2.3 and 2.7 show substantially greater levels of KH-induced TKE aloft than those observed in the lower stability, lee-wave response scenarios. At *x* = −3 km, both strong stability cases show substantial first-layer stagnation zone turbulence and some level of secondary IGW breaking above. At *x* = 0 km, the maximum TKE in the stagnation zone is higher than that of the turbulent ABL below. Farther downstream at *x* = 2 and 4 km, the stagnation zone and turbulent ABL have combined and subsumed the shooting flow in a hydraulic-jump-like response.

Figure 6 presents shear generation and buoyant production of turbulence for the Fr^{−1} = 1.6 and 2.3 cases of the trapped-lee-wave and nonlinear-breakdown regimes, respectively. Close inspection of these plots shows the onset of stagnation zone turbulence is due to both shear generation (between the upper edge of the shooting flow and the lower edge of the stagnation zone) and buoyant production in the upper portion of the stagnation zone, with each being an aspect of the stagnation zone KH instability. In contrast to the stagnation zone, the leeside boundary layer TKE growth is mainly through shear production in these cases where the turbulent ABL is near-neutrally stratified; however, for the leeside nonlinear-breakdown case, buoyant sink–source is observed downstream, where mixing of air from above creates local pockets of convective stability–instability.

These two cases (similar regime cases omitted for clarity) are representative of the general trend that, as stability (or potential energy of the scenario) is increased, there is enhanced buoyant production in the stagnation zone. In addition, as stability increases, shooting flow leads to greater shear production at stagnation zone upper and lower height bounds. Since the wavelength of the vertically propagating MW decreases as Fr^{−1} increases, the vertical velocity gradient in the wind profile increases, leading to the larger levels of shear production evident in Fig. 6. This trend in growing MW activity associated with vertical wavelength is consistent with previous MW efforts as far back as Queney (1948). These LES results show that, when the vertical wavelength becomes comparable to the sufficiently energetic turbulent scales in the stagnation zone and ABL, an abrupt change in leeside flow response occurs, resulting in transition to the hydraulic-jump-like state.

The shear and buoyancy influences shown in Fig. 6 show considerable agreement with the results of Smith and Skyllingstad (2009) in that buoyancy perturbations manifest primarily as a source of turbulence in the stagnation zone but act as a sink at the stagnation zone–shooting-flow interface and arise as both source and sink contribution locally within downstream rotors. In addition, our results show agreement with respect to the substantial shear generation of turbulence at the stagnation zone–shooting-flow interface and within rotor structures. An important difference to note is that the neutrally stratified ABL here indicates these effects are not due strictly to heating or cooling at the surface, but rather to the mere presence of a fully turbulent ABL. Nonetheless, we fully expect surface heating to influence ABL stability and turbulence; thus, further investigation as to the combined impact on MW development and leeside flow response regime transition would be beneficial.

## 4. Compensated Froude characterization and the importance of a turbulent ABL

Classical linear theory has been extensively used to characterize flow behavior regimes for flow splitting (in the case of a three-dimensional hill) and IGW breaking based on Fr^{−1} (Smith 1989). The regime diagram of Smith (1989) provides the linear estimates for Fr^{−1} at which the transition from linear wave into the nonlinear wave (IGW breaking) regime will occur for a given horizontal mountain aspect ratio *r* = *a*_{y}/*a*_{x}, where *a*_{y} and *a*_{x} are the spanwise and streamwise hill widths, respectively. For the two-dimensional flows analyzed here, we consider the asymptotic nature of the linear theory estimates for *r* → ∞. In this case, the onset of the IGW-breaking regime as estimated by linear theory corresponds to Fr^{−1} ≈ 0.6, from Fig. 5 in Smith (1989). However, the author provides in the regime diagram values from nonlinear calculations suggesting an increase in the transition value by 30% such that the nonlinear value is Fr^{−1} ≈ 0.78. It is worth noting that the linear theory may not be applicable near the transition, where the emergence of stagnation (nonlinearities) in the flow would invalidate its estimates. In addition, neither the linear theory estimates nor the nonlinear-calculated values account for the presence of a turbulent, and by definition nonlinear, atmospheric boundary layer. From the results shown in the previous sections, the presence of a turbulent ABL has a substantial impact on the dynamics of MW flows. Therefore, we propose a compensation for this omission in classical Froude formulation through a more representative choice of vertical length scale in the Froude number definition. Taking the difference between hill height and boundary layer height for the characteristic vertical length scale, we define an inverse compensated Froude number as

This choice of vertical length scale is motivated by the influence of the boundary layer in reducing the effective terrain height as perceived by the stably stratified flow above. When *h* ≤ *z*_{i}, the compensated Froude number becomes zero or negative, indicating the ABL presence effectively hides the hill from the upper-level stratified flow. This is in general agreement with the notion that surface friction processes attenuate MW activity, as suggested by Richard et al. (1989), Ólafsson and Bougeault (1997), and Smith (2007). The corresponding value of inverse compensated Froude number for the transition between linear MWs and nonlinear, IGW-breaking MWs suggested by the simulations presented in the previous sections, , differs substantially from the suggested nonlinear-calculated transition value of 0.78 from Smith (1989). As further evidence of the advantage in choosing a vertical length scale accounting for turbulent boundary layer height, Fig. 7 shows instantaneous zonal velocity and TKE contours (as in Fig. 3 and Fig. 4) for three scenarios with identical Fr^{−1} = 1.2, but different under modified hill height and Brunt–Väisälä frequency with constant boundary layer height (*z*_{i} = 500 m). The = 0 scenario shows a hidden hill (i.e., *h* = 315 m < *z*_{i}) with very strong stratification (*N* = 0.02 s^{−1}) and no IGW breaking. The = 0.2 case has a moderate hill height (*h* = 600 m) just taller than the boundary layer and under moderate stability aloft (*N* = 0.02 s^{−1}) and again shows no IGW breaking. The bottom row has = 0.7, corresponding to a tall hill (*h* = 1200 m) under the weakest simulated stability (*N* = 0.01 s^{−1}), and exhibits IGW breaking. For all of these scenarios, the classical Fr^{−1} = 1.2 would imply IGW breaking; however, the compensated Froude appropriately predicts a dearth of IGW breaking in the two scenarios where < 0.5 and the presence of IGW breaking for > 0.5. Additionally, all of the simulations investigated in this manuscript (including those of the previous sections and Fig. 7) suggest ≈ 1.0 is the appropriate regime transition value between trapped-lee-wave and nonlinear-breakdown states. We fully expect that stability of the boundary layer itself would increase/decrease the effective for stable–convective ABLs. In fact, considering the results demonstrating this aspect of MW flows in Smith and Skyllingstad (2009), an additional extension to the compensated Froude number proposed herein should be included to account for this sensitivity.

One important aspect of atmospheric LES investigations is the premise that resolved turbulence must exist in both boundary and initial conditions. Figure 8 presents streamwise velocity and TKE of two simulations (Fr^{−1} = 1.6 and 2.3) in which inflow boundary and initial conditions for streamwise velocity were constant with height and void of resolved turbulence. Several key differences between these simulations and the corresponding simulations with a turbulent ABL (see Figs. 3 and 4) are found: (i) an upstream ABL interacts with upstream upper-level flow, producing variability in *z*_{i} and subsequently horizontal variability aloft, instigating more TKE in the stagnation zone; (ii) more TKE generated in the stagnation zone results in a more prominent critical layer that constrains the upward propagation of trapped lee waves and increases dissipation from above; (iii) the initial presence of a downstream, leeside ABL provides resistance to the development of the shooting flow, results in smaller and weaker rotors, and increases dissipation from below; (iv) persistent and fully developed turbulence is achieved in both the stagnation zone and ABL regions only when resolved turbulent structures are included in both boundary and initial conditions. From these results, it is clear that any atmospheric LES investigation of MW development must include a fully developed three-dimensional turbulent boundary layer in order to realistically capture the complex interactions between stability-induced MW flow features and the ABL. Earlier studies using mesoscale models with one-dimensional, fully parameterized turbulence models found inclusion of an ABL suppresses gravity wave breaking aloft (see Smith 2007). Recent results by Muñoz-Esparza et al. (2016) indicate this is strictly due to the misrepresentation of three-dimensional turbulence effects by one-dimensional parameterization schemes. Also, although not explicitly mentioned in Doyle et al. (2011), the authors include comparison of an LES simulation without a resolved, turbulent-inflow ABL against a free-slip mesoscale simulation, wherein nearly identical results were obtained between the two modeling approaches.

## 5. Energy spectra and turbulent structures

Turbulent kinetic energy spectra provide a quantitative analysis of energy content across a range of wavelengths. To further elucidate the characteristics of disparate leeside flow response regimes, TKE spectra at two locations of similar relevance across scenarios are presented (see red markers in Fig. 4). The left panel of Fig. 9 shows the TKE spectra at the height of minimum streamwise velocity in the stagnation zone [], directly above the downstream shooting-flow maximum on the lee side of the hill (which is Froude number dependent). The right panel shows the TKE spectra at a height , with the downstream location corresponding to the secondary maximum of streamwise velocity (again Froude number dependent). These two locations are considered relevant here, because in the former case stagnation zone turbulence is established, and in the latter the location is where the interaction between shooting flow and the turbulent boundary layer occurs. The coordinates of these locations vary depending on to account for decreasing characteristic vertical wavelength and the reduction of leeside streamwise wavelength with increasing inverse Froude. The inflow ABL TKE spectra at *z* = *z*_{i}/2 is plotted in both panels as a reference for the equilibrium inflow boundary layer turbulence.

For = 0.4, as shown in previous results, development of TKE in the stagnation zone is negligible. However, at the downstream boundary layer location, the TKE spectra show the enhancement of TKE above the equilibrium reference resulting from the descending weak shooting flow and subsequent increase in shear production of TKE. For this weakest stability case, the downstream ABL spectrum clearly shows a near-uniform shift in energy content across all scales relative to the equilibrium ABL (i.e., the curvature of the spectra with wavenumber is nearly identical). As increases into the trapped-wave-response cases (0.55, 0.7, and 0.8), turbulence in the stagnation zone increases substantially. Interestingly, energy spectra over the three trapped-lee-wave cases are extremely similar in energy content and distribution. This pattern holds true for the boundary layer spectra as well. For the leeside nonlinear-breakdown regime cases, = 1.0, 1.2, and 1.5, both the stagnation zone and downstream ABL spectra are significantly more energetic across all scales than the preceding lee-wave energy spectra. Despite spanning a range of 0.5 in , these spectra again nearly overlap. It is important to note again that the locations of the points in each case are different; however, the intent was to choose points that, although different in coordinates, are phenomenologically similar across and representative of particular interactions in the flow (i.e., maximum stagnation zone turbulence production and strong interaction between shooting flow and boundary layer). The similarity of energy spectra between scenarios of similar flow type supports the existence of an abrupt transition in regime parameter space and suggests two distinct energy distribution states. It appears that flow response is closely tied to a particular relationship between stagnation zone and ABL turbulence levels, leading to highly similar spectral distributions between Froude number cases of a given regime. In an attempt to further constrain the threshold leeside response transition value, a simulation for = 0.9 was performed (not shown), which resulted in the leeside flow response of a strongly dissipated trapped lee wave. From this additional result, along with the results shown in previous sections, the transition value for leeside flow regime is in the range 0.9 < < 1.0.

As a qualitative, yet poignant illustration of the complex role turbulence plays in both the stagnation zone and ABL, the interactions of turbulence in these locations, and the subsequent impact on MW flow regimes, Fig. 10 shows an isovolume rendering of enstrophy *ξ*, or the square of vorticity magnitude (see Kundu and Cohen 2004), for ξ = 0.02 s^{−2}, colored by potential temperature for = 0.4, 0.8, 1.0, and 1.5. In this figure, it is again observed that the weakest stability case produces no stagnation turbulence aloft, but an enhancement of boundary layer turbulence occurs owing to the increased shear from the weak shooting flow amplifying shear generation in the boundary layer. For the trapped-lee-wave scenario of = 0.8, moderate turbulence levels and vertical wavelength allow the essentially laminar shooting flow to maintain a clear separation between upper-level critical layer and lower-level boundary layer turbulence. This can be seen from the lack of mixing between warm turbulent air aloft and the cooler turbulent ABL air near the surface until farther downstream, where the lee wave is dissipated by the turbulence in the two regions. The empty space between is precisely where the nearly laminar (zero enstrophy) shooting flow remains trapped by the two turbulent layers. In the two leeside nonlinear-breakdown cases, = 1.0 and 1.5, mixing of the two layers occurs with increasing efficiency at shorter downstream distance as inverse Froude grows. This is easily seen from the weak warm air mixing into the near surface in the = 1.0 case and substantial mixing of warm air into the near surface in the = 1.5 case. Also, progressive enhancement of second-level wave breaking with increasing is observed, with a clear separation from the first-level stagnation zone and ABL turbulence regions. These qualitative visualizations support the notion that nonlinear dynamics play a critical role in the complex interactions between stably stratified flows over terrain and a turbulent ABL, requiring numerical models that allow explicit resolution of such phenomena in order to properly model these types of complex scenarios.

## 6. Conclusions

The focus of this work is to provide new insights into the dynamics of a turbulent ABL on MW development and regime transitions, namely the onset of IGW breaking and leeside flow response. To that end, we perform simulations of large-scale stably stratified flow over an isolated, two-dimensional hill using turbulence-resolving large-eddy simulation. Our results demonstrate that a turbulent ABL influences MW development in critical aspects, such as dissipation of trapped lee waves and amplified turbulence. In particular, enhanced stagnation zone turbulence is triggered by both upstream boundary layer undulations interacting with upper-level air and the spatiotemporal heterogeneity of downstream interactions between the shooting flow and turbulent ABL. The heterogeneity of these interactions propagates vertically and upstream, exhibiting a positive feedback cycle on the generation of three-dimensional turbulence through KH instability in the stagnation zone.

An inverse compensated Froude number, = *N*(*h* − *z*_{i})/*U*_{g}, with vertical length scale taken as the difference between maximum hill height and boundary layer height, is proposed to account for the influence of the turbulent ABL on the transition between both nonbreaking and breaking IGW, as well as the leeside flow response regime from no response to trapped-lee-wave to nonlinear-breakdown (or hydraulic-jump-like flow) states. However, while our compensated Froude number accounts for upper-level stratification influence through the Brunt–Väisälä frequency, it accounts for neither capping inversions, nor vertical variability in stratification aloft (*N* is constant here). Further investigation is required to systematically evaluate such dependencies, the former of which has been shown from observations (Mobbs et al. 2005) and LES modeling (Smith and Skyllingstad 2011) to have an important role in determining the formation and strength of MW activity. Also, considering the results in Smith and Skyllingstad (2009), we fully expect that stability of the boundary layer will itself influence regime transition, and the compensated Froude number proposed herein should be extended in order to account for this sensitivity in the definition of MW regimes. Additionally, and although not presented here, exploratory simulations performed with different nondimensional mountain half-widths showed influence on the onset of IGW breaking and leeside flow response, making this yet another critical aspect to be considered when proposing a more general Froude number definition.

Through analysis of vertical profiles of TKE and budget terms for shear production and buoyancy source–sink on the lee side of the hill, the interactions leading to various flow regimes are elucidated in detail. First, the onset of trapped lee waves appears when the vertical wavelength provides separation between the KH-induced stagnation zone (critical layer) and the boundary layer turbulence so as to constrain the oscillatory ascension and descension of the leeside shooting flow. Enhanced dissipation occurs through lee-wave interactions at lower and upper interfaces with one or both of these turbulent regions. Over a range of inverse compensated Froude numbers from 0.55 to 0.9, trapped lee waves develop and undergo increasing levels of dissipation with downstream distance as inverse Froude increases. Transition to the nonlinear-breakdown regime is shown here to occur for ≈ 1.0. This transition is driven by a combination of enhanced turbulence in both the stagnation zone and turbulent boundary layer in addition to the decreasing vertical wavelength of IGW.

These conclusions should be regarded with care, given the idealization of our environmental forcing. This forcing has been recently found by Hills et al. (2015) using a mesoscale-like model to be a key factor in determining the preference for trapped lee waves to either leak into the stratosphere or dissipate into the boundary layer. Nevertheless, as recently demonstrated by Muñoz-Esparza et al. (2016), mesoscale models using one-dimensional PBL parameterizations result in unphysically persistent, trapped lee waves, owing to a misrepresentation of three-dimensional turbulence effects and interaction with the large-scale flow. This is once again another example where LES modeling can help in providing further understanding of the fundamental aspects of complex terrain–atmosphere interactions.

Turbulent kinetic energy spectra were analyzed at phenomenologically similar locations across the range of Froude scenarios investigated here. These spectra exhibited remarkably similar distribution of energy content across all scales for Froude cases of a given regime. In the trapped-lee-wave cases of = 0.7 and 0.8, the energy spectra in the stagnation zone and ABL locations of interest were nearly indistinguishable across Froude values. For a larger range of stability conditions, the leeside nonlinear-breakdown cases of = 1.0, 1.2, and 1.5 showed overlapping spectra in the stagnation zone as well as at the downstream ABL location, where interactions between ABL and shooting flow were expected to be strongest. This analysis of energy spectra supports and extends the understanding of the abrupt transition between flow response regimes at some critical value of inverse compensated Froude number of 0.9 < < 1.0. Finally, the complex role of turbulence in both the stagnation zone and boundary layer inherent to different MW regimes was conveyed through qualitative visualization of enstrophy isovolumes. These results support the notion that nonlinear dynamics play a critical role in the complex interactions between stably stratified flows over terrain and a turbulent ABL.

While the efforts here advance the understanding and state of the art in LES modeling of MW flows, particularly when the near-surface region is the targeted area of interest, there is more to be explored. Extension of these efforts to isolated three-dimensional hills should be performed in order to gauge the effect of a turbulent ABL on flow splitting. Furthermore, the effects of Coriolis-induced veering should also be reconsidered in the three-dimensional and explicitly resolved turbulent ABL setting. Finally, and potentially most critically, the efforts here provide ample motivation and demonstrated need for substantial improvement in the spatial and temporal fidelity of observational campaigns in order to facilitate validation and comparison with sophisticated LES modeling results in the context of real-world complex terrain and forcing conditions.

## Acknowledgments

The majority of efforts presented here and in preparation of this manuscript were carried out during the postdoc tenures of JAS and DME in the Earth and Environmental Sciences Division at Los Alamos National Laboratory. This research was supported by the Laboratory Directed Research and Development (LDRD) program at Los Alamos National Laboratory (20130487ER). This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract DE-AC52-06NA25396. The authors thank the three anonymous reviewers for their thorough and constructive suggestions to improve the manuscript. We would also like to thank Dr. François Pimont for his insightful advice regarding the large-scale pressure gradient force parameterization used in this work.

## REFERENCES

*Proc. IBM Scientific Computing Symp. on Environmental Sciences*, White Plains, NY, International Business Machines Corporation, 195–210.

*Bull. Amer. Meteor. Soc.*,

**34**, 205–211.