The internal structure and dynamics of rotors that form in the lee of topographic ridges are explored using a series of high-resolution eddy-resolving numerical simulations. Surface friction generates a sheet of horizontal vorticity along the lee slope that is lifted aloft by the mountain lee wave at the boundary layer separation point. Parallel-shear instability breaks this vortex sheet into small intense vortices or subrotors.
The strength and evolution of the subrotors and the internal structure of the main large-scale rotor are substantially different in 2D and 3D simulations. In 2D, the subrotors are less intense and are ultimately entrained into the larger-scale rotor circulation, where they dissipate and contribute their vorticity toward the maintenance of the main rotor. In 3D, even for flow over a uniform infinitely long barrier, the subrotors are more intense, and primarily are simply swept downstream past the main rotor along the interface between that rotor and the surrounding lee wave. The average vorticity within the interior of the main rotor is much weaker and the flow is more chaotic.
When an isolated peak is added to a 3D ridge, systematic along-ridge velocity perturbations create regions of preferential vortex stretching at the leading edge of the rotor. Subrotors passing through such regions are intensified by stretching and may develop values of the ridge-parallel vorticity component well in excess of those in the parent, shear-generated vortex sheet. Because of their intensity, such subrotor circulations likely pose the greatest hazard to aviation.
Stratified airflow over mountainous terrain can lead to a rich spectrum of atmospheric responses over scales ranging from planetary to turbulence. One of the most severe topographically forced phenomena known is the rotor, which is a low-level vortex with a circulation axis oriented parallel to the mountain ridgeline. The surface winds associated with rotors are typically characterized by strong downslope winds near the surface that decelerate rapidly in the lee and give way to a weaker recirculating flow directed back toward the mountain (e.g., Homboe and Klieforth 1957; Kuettner 1959; Doyle and Durran 2002). Rotors pose a substantial aeronautical hazard owing to the potential for very strong lower-tropospheric turbulence and shear, and have been suggested to contribute to numerous aviation incidents and accidents (e.g., NTSB 1992; Lester 1993; Carney et al. 1996; Darby and Poulos 2006). For example, a severe turbulence incident, likely associated with a low-level rotor, resulted in the loss of an engine on a commercial United Airline Boeing 747–100 at 600-m AGL near Anchorage, Alaska (Kahn et al. 1997). In spite of their clear significance to the meteorology and aviation communities, the dynamics and structure of rotors are poorly understood and forecasted, in part because of infrequent and insufficient observational measurements, and inadequate sophistication and fidelity of numerical weather prediction models.
Mountain waves and rotors were the subject of two of the first modern U.S. multiagency field programs in meteorology, the Sierra Wave Project (SWP) and its follow-on, the Jet Stream Project (JSP), both of which took place in the early 1950s (Holmboe and Klieforth 1957; Grubišić and Lewis 2004). With the exception of research aircraft observations of several rotor events in the lee of the Rocky Mountains (Lester and Fingerhut 1974) and occasional serendipitous remote sensing lidar measurements of rotors (Banta et al. 1990; Ralph et al. 1997), there were remarkably few direct observations of rotors during the first four decades following the SWP and JSP. This situation has now changed as rotors have become the focus of new observational campaigns.
In a recent study, the near-surface flow across and downwind of the Wickham mountain range on the Falkland Islands was observed during a field campaign aimed at improving the prediction of orographically generated turbulence (Mobbs et al. 2005). Several strong downslope windstorm events, punctuated by episodes of short-lived periods of flow separation and rotor development, were documented in the lee of the ridge crest. Darby and Poulos (2006) documented the evolution of a mountain wave/rotor system interacting with an approaching cold front in the lee of Pike’s Peak in the Rocky Mountains using Doppler lidar, wind profiler, and research aircraft observations, as well as numerical simulations. They found that the wind shear associated with mountain waves and rotors evolved rapidly as a result of the mesoscale and synoptic conditions modulating the upstream flow properties. The Sierra Rotors Project (SRP), which took place in the spring of 2004, consisted of a suite of ground-based observing platforms including a network of surface stations, two wind profilers, and radiosonde observations upstream and downstream of the Sierra Nevada range. The SRP measurements reveal a number of events with accelerated downslope flow along the lee slope together with reversed flow and rotors further downstream (Grubišić and Billings 2007). Downslope windstorm events during the SRP were observed to occur most frequently in the local afternoons due to diurnal boundary layer heating effects that are manifested as a multiscale dynamical response (Jiang and Doyle 2005). In the Terrain-induced Rotor Experiment (T-REX), which took place in the Sierra Nevada range in the spring of 2006, the structure and evolution of atmospheric rotors were explored using a comprehensive suite of ground based, remote sensing, and airborne platforms (Grubišić et al. 2004). Unfortunately, none of the above field programs, with the exception of the just-completed T-REX program, have provided sufficient systematic and detailed measurements of the internal structure of rotors to establish the nature of smaller-scale circulations and turbulence within rotors, which is the primary focus of this study.
Numerical simulations, conducted by Doyle and Durran (2002), have suggested that a key aspect of rotor development involves the mutual interaction between the lee wave and the surface boundary layer. Their two-dimensional simulations indicate that a thin vortex sheet, generated by mechanical shear in the boundary layer, separates from the surface owing to adverse pressure gradients associated with lee waves. To explore the dependence of rotors on the lee wave amplitude, Doyle and Durran (2002) also conducted a series of simulations with varying mountain heights and interface depths in a background flow with a two-layer stratification. The results indicate that the magnitude of the reversed flow in the primary rotor for a simulation with surface friction is highly correlated with the strength of the adverse pressure gradient in the lee wave in an otherwise identical simulation without surface friction (i.e., with a free-slip lower boundary condition). Other studies have also confirmed the link between the lee wave amplitude and rotor characteristics. Vosper (2004) conducted a series of simulations and found that, as the ratio of the mountain height to inversion height (H/zi) increases, the lee wave amplifies and exceeds a critical threshold for the onset of flow separation triggering the development of a rotor. Vosper found that the flow state transitioned between lee waves, rotors, and hydraulic jumps as a result of changes to a sharp temperature inversion located near the mountain crest level, which impacts the upstream Froude number, Fi, for a two-layer stratification. Similarly, Mobbs et al. (2005) identified Fi and H/zi to be the two key parameters delineating the windstorm and rotor regimes for the Falkland Islands field campaign. The two-dimensional simulations of Hertenstein and Kuettner (2005) underscore the significance of wind shear within low-level inversions. Hertenstein and Kuettner identified two distinct flow states associated with the characteristics of the vertical shear within the inversion that they refer to as type 1 and type 2 rotors. A recirculating rotor forms beneath the lee wave crest in the presence of forward shear in the low-level inversion. When the shear within the inversion layer is reduced, Hertenstein and Kuettner hypothesized that a type 2 rotor forms, which has similarities to an unsteady wave breaking state (e.g., Afanasyev and Peltier 1998) or a hydraulic jump (Lester and Fingerhut 1974; Dawson and Marwitz 1982).
The numerical simulations of Doyle and Durran (2002) and Hertenstein and Kuettner (2005) both suggest that vertical and horizontal wind shear and turbulence production are maximized along the elevated sheet of horizontal vorticity, particularly along the upstream side of the lee wave or leading edge of the rotor circulation, in general agreement with anecdotal evidence from the SWP (Holmboe and Klieforth 1957) and the aircraft observations of rotors taken in the lee of the Rocky Mountain Front Range analyzed by Lester and Fingerhut (1974). During the JSP, an instrumented B29 aircraft penetrated a rotor circulation and encountered severe turbulence. The flight-level data shown in Fig. 1 indicates that the aircraft encountered a number of strong vertical gusts approaching ±20 m s−1 over a 50-s period. In another event during the JSP, an instrumented glider was destroyed in midflight after an encounter with a rotor (Holmboe and Klieforth 1957). These rare encounters with rotors suggest the presence of extreme turbulence that is likely composed of intense smaller-scale vortices, referred to here as “subrotor” circulations. Although no comprehensive observations of subrotors exist, scanning lidar observations of low-level airflow during the 9 January 1989 severe downslope windstorm over Boulder, Colorado, exhibit smaller-scale circulations resembling subrotor eddies, with a horizontal scale of ∼1 km, that are embedded within a recirculating rotor flow positioned beneath a lee wave (Banta et al. 1990). Time lapse photographs of rotors (e.g., Ozawa et al. 1998) are also suggestive of smaller-scale circulations embedded within topographically forced rotors.
Although it has been suggested that the most severe turbulence occurs along the leading edge of rotors, considerable uncertainty exists with regard to the structure and nature of the turbulence within rotors. The internal structure of the rotor remains relatively unexplored, in part because of the hazardous nature of the turbulence that has prevented systematic in situ observations and in part because the computational power to explicitly resolve the eddy structure within the rotor has only recently become available. Rotor simulations that make use of sufficient horizontal resolution (Δx ∼ 100 m) to capture the internal rotor structure have been generally limited to two-dimensional models (e.g., Doyle and Durran 2002; Hertenstein and Kuettner 2005). Previous three-dimensional simulations of rotors generally have used horizontal grid increments of 300 m or greater, which cannot adequately resolve the eddy structure within the rotor (Doyle and Durran 2004; Grubišić and Billings 2007; Darby and Poulos 2006).
Although they are computationally efficient, 2D rotor simulations may yield misleading results. Investigations of constant wind speed and static stability flow past elongated 3D mountains suggest that the response can differ significantly relative to flow over 2D ridges (Epifanio and Durran 2001). Flow splitting, vortex shedding, and wakes are other potentially relevant phenomena that may occur in stratified flows past 3D obstacles, but are not captured in a 2D geometry (e.g., Smolarkiewicz and Rotunno 1989; Schär and Durran 1997; Epifanio and Durran 2002a, b). The sensitivity of rotor characteristics to three-dimensional aspects of the topography has yet to be established.
In this study, very high resolution numerical simulations of rotors are performed. The primary objectives of this study are to (i) contrast the characteristics of rotors in 2D and 3D flows, (ii) explore the impact of three-dimensional terrain variations on rotor characteristics and their internal structure, and (iii) identify the sources of vorticity in rotors and subrotors. In the following section the numerical model is described. The results of 2D simulations are presented in section 3. Section 4 contains a discussion of 3D simulation results. The instability mechanism for small-scale circulations within the rotor is described in section 5. The summary and conclusions are presented in section 6.
2. Numerical model description
The atmospheric portion of the Naval Research Laboratory’s Coupled Ocean–Atmospheric Mesoscale Prediction System (COAMPS; Hodur 1997) is used to conduct these simulations. A brief overview of the COAMPS model is presented here. The prognostic variables include the Cartesian wind components (ui, where i = 1, 2, 3), perturbation Exner function (π′), and potential temperature (θ); the effects of moisture are neglected. The prognostic equations expressed using tensor notation are
where cp is the specific heat at constant pressure, Rd is the dry gas constant, g is the acceleration due to gravity, ρ is the density, c is the speed of sound, ν is the hyperdiffusion coefficient, and f the Coriolis force. In this study, the Coriolis force is specified as 10−4 s−1. The overbar variables correspond to the reference state, which is horizontally homogeneous and in hydrostatic balance. The turbulent subgrid-scale fluxes for momentum, τij, and heat, Hj are included. Terrain is incorporated through a transformation to the following coordinate:
where zt is the depth of the model computation domain, z is the physical height, h is the terrain elevation, and σ is the transformed vertical coordinate.
Equations (1)–(3) are discretized on an Arakawa C grid using a Cartesian projection. The equations are efficiently integrated using a time splitting technique that features a semi-implicit treatment for the vertically propagating acoustic waves following Klemp and Wilhelmson (1978). The time differencing is centered for the large time step. The horizontal advection and horizontal smoothing terms are represented by fourth-order accurate differencing, and second-order differencing is used to represent the vertical advection, pressure gradient, and divergence terms. A Robert time filter is applied to all predicted variables.
The subgrid-scale turbulent flux parameterization is based on Lilly (1962) and expressed as
where Dij is the deformation tensor and KM and KH are the eddy viscosity and diffusion coefficients computed using the Lilly–Smagorinsky local equilibrium model (Lilly 1962; Smagorinsky 1963). The eddy diffusivity for momentum and heat are specified as
where λ = csΔ, cs = [c3m/(cɛ1 + cɛ2)]1/4, and the mixing length, Δ, is
The deformation, S2, is defined as
The values of the coefficients are based on Stevens et al. (1999) such that cm = 0.0856, ch1 = cm, ch2 = 0.1184, cɛ1 = 0.1911, and cɛ2 = 0.6539. The general implementation of the subgrid-scale turbulence parameterization is described in Golaz et al. (2005). The surface momentum flux is computed following the Louis (1979) and Louis et al. (1982) formulation, which makes use of Monin–Obukhov similarity theory and assumes a vegetation roughness length of 1 cm. The surface heat flux is zero, implying that the surface ground temperature is assumed to be in balance with the surface air temperature.
The radiation condition proposed by Orlanski (1976) is used for the lateral boundaries with the exception that the Doppler-shifted phase speed (u ± c) is specified and temporally invariant at each boundary (Pearson 1974; Durran et al. 1993). The mitigation of reflected waves from the upper boundary is accomplished through a radiation condition formulated following the Durran (1999) approximation to the Klemp and Durran (1983) and Bougeault (1983) techniques.
The 2D and 3D simulations are initialized using a reference state approximating the conditions upstream of the Colorado Front Range on 3 March 1991. This particular case is significant because, not only were rotors observed, but also because a United Airlines B737 crashed while attempting to land at Colorado Springs during the event. The sounding used to initialize the basic state in these simulations, as shown in Fig. 2, is based on the potential temperature and cross-mountain wind component from the Grand Junction and Denver, Colorado, and Lander, Wyoming, radiosondes at 1200 UTC 3 March 1991. The cross-mountain wind component throughout the troposphere increases from roughly 10 m s−1 at the surface to 45 m s−1 near the top of the domain. The potential temperature profile comprises a statically stable layer near the surface with N ∼ 0.009 s−1, capped by an inversion in the 2500–3300-m layer with N ∼ 0.017 s−1. The model was initialized with small random and spatially uncorrelated potential temperature perturbations of maximum magnitude 0.1 K that were introduced in approximately the lowest 3 km of the model.
The topography for the 2D simulations, and for 3D simulations using an infinitely long 2D ridge, are specified using a two-sided Witch of Agnesi profile:
with a = au for x < 0 and a = ad otherwise. The Front Range is asymmetric with the lee slopes steeper than the upwind slopes. To represent this asymmetry, the upwind half-width is set to au = 15 km and the downwind half-width to ad = 5 km. The mountain height is hm = 1500 m, which is approximately representative of the mean elevation gain between the eastern plains and the crest of the Colorado Front Range. The surface pressure far away from the ridge axis is specified to be equivalent to the standard atmospheric pressure at an elevation of 1.5 km MSL.
The fully 3D topography is a smooth ridge following Epifanio and Durran (2001). The long axis is oriented normal to the x direction and is defined
The horizontal aspect ratio (the ratio of the y to x dimension of the ridge) is β = 4, and hm = 1500 m. Consistent with the two-dimensional simulations, the upwind (downwind) half-width is 15 km (5 km). In several numerical experiments an isolated peak protrudes above the ridgeline as illustrated in Fig. 3. This peak is obtained by adding a second mountain of the form (12) to the topography with hm = 500 m, β = 1, au = 7.5 km, and ad = 3.25 km in which the origin is shifted “north” by 20 km.
To achieve sufficient resolution to explicitly resolve the internal structure of the rotors, horizontally nested grids are utilized in the 2D and 3D simulations. The 2D simulations use three grid meshes with horizontal resolutions of 540, 180, and 60 m, with 601, 241, and 397 points on each mesh, respectively. The 3D simulations using an infinite 2D ridge (in the y direction) have three grid meshes with horizontal resolutions of 540, 180, and 60 m corresponding to horizontal dimensions of 601 × 25, 241 × 55, and 397 × 139 points on each mesh, respectively. The 3D simulations using a finite 3D ridge have four grid meshes with horizontal resolutions of 1620, 540, 180, and 60 m, corresponding to horizontal dimensions of 223 × 223, 205 × 343, 145 × 199, and 217 × 217 points on each nested grid mesh, respectively. The position of the 3D topography on the 540-m mesh is shown in Fig. 3. In both 2D and 3D simulations a total of 90 vertical levels are used, and the vertical grid spacing stretches from 15 m at the lowest level to a constant 50 m in the layer between z = 250 and 2550 m and then gradually stretches to a maximum spacing of 400 m at the model top near 12.3 km.
3. Infinitely long uniform ridge
In this section we compare rotors generated by a purely 2D flow over a ridge whose cross section is defined by (11) with those in a closely related 3D geometry. In the 3D case, the domain is periodic in the spanwise (y) direction, and the mountain profile and upstream atmospheric conditions are independent of y and identical to those in the purely 2D case. The large-scale wave and rotor structure for the purely 2D simulation is shown in Fig. 4, which is a plot of isentropes and horizontal wind speed for a 135-km-wide subdomain of the coarsest (540 m) resolution mesh at t = 3.5 h. A series of trapped mountain waves are apparent in the lee of the terrain. The waves are trapped by the increase of the cross-mountain winds with height and by the decrease in static stability above the inversion layer, apparent near 3 km in Fig. 2 (Scorer 1949). The peak-to-trough displacement of the isentropes in the lee waves is approximately 1.8 km. A region of reversed recirculating flow is present beneath each lee wave crest. The series of rotor circulations beneath each crest is qualitatively similar to the rotor streaming flow discussed by Förchtgott (1949).
Time-mean fields from the highest resolution (60 m) nest, averaged over the interval between 3 and 3.75 h, are shown in Fig. 5. The top panels show fields from the purely 2D case, whereas the bottom panels show fields from the 3D case that have been averaged over the 8280-m “north–south” extent of the 60-m nested mesh. In both simulations, the winds, isentropes, and the distribution of y-component vorticity (η) are very similar outside and along the upstream edge of the rotor, but there are marked differences within the interior of the rotor itself. In the 2D simulation, the speed of the time-averaged near-surface reversed flow exceeds 26 m s−1, whereas it is only about 7 m s−1 in the 3D case. The sheet of intense η generated by boundary layer friction that lifts off the surface at the leading edge of the lee wave/rotor also evolves very differently in the 2D and 3D cases. In 2D the vorticity is only subject to weak dissipation and the time-averaged flow within the rotor contains a large region of positive η. A smaller region of negative η is also present, again generated by boundary layer friction, but in this case by the strong reversed flow at the base of the rotor. On the other hand, in the 3D case, the sheet of vorticity rising along the upstream edge of the lee wave is dissipated more rapidly so that the y-component vorticity of the time- and y-averaged velocity field within the rotor is very weak. The contrast between the vorticity structures inside the main rotor in Figs. 5b and 5d is a consequence of the well-known fundamental difference between 2D and 3D turbulence, namely that there is a systematic cascade of energy to small scales in 3D, but not in 2D (Tennekes 1978).
The evolution of the y component of the horizontal vorticity for the 2D case is shown in Fig. 6 by three snapshots taken at 1-min intervals beginning at 3.5 h using data from the highest resolution (60 m) mesh. In contrast to the time-averaged picture in Fig. 5b, a much more complex structure is apparent. As the vortex sheet is lifted off the surface by the lee wave, it breaks up into several individual vortices, or subrotors. The instability that causes the vortex sheet to break up into to these subrotors will be explored further in section 5; here we will focus on the behavior of the subrotors after their generation. Three individual subrotors are identified as s1, s2, and s3 in Fig. 6. Subrotor s2 has just formed, due to the rollup of the main vortex sheet in Fig. 6a, and then is subsequently advected upward along the leading edge of the lee wave/large-scale rotor in Figs. 6b and 6c. Subrotor s1 begins as just a local maximum of shear vorticity (Fig. 6b) but subsequently intensifies as its section of the vortex sheet begins to role up (Fig. 6c). Subrotors such as s1 and s2 appear to form along the vortex sheet at a frequency of approximately every 2–3 min. Subrotor s3 has already advected over the crest of the lee wave at the time shown in Fig. 6a and continues to descend and become entrained into the large-scale rotor circulation in Figs. 6b and 6c. The vorticity associated with the circulation inside the large-scale rotor appears to be continually replenished by the gradual dissipation of the individual subrotors.
The preceding 2D results may be compared with the evolution of the y component of the horizontal vorticity in the 3D case shown in Fig. 7, which shows contours of η on the 60-m resolution mesh in a vertical section normal to the ridge along the centerline of the grid at 1-min intervals beginning at 3 h 30 m. As in the 2D case, a series of subrotors develop along the vortex sheet and advect downstream. However, in contrast to the 2D simulation, there is no systematic transfer of positive η vorticity between the subrotors and the interior of the main large-scale rotor. Instead, most of the subrotors are swept downstream in the flow along the top of the main lee wave/rotor, and the interior of the large-scale rotor is filled with almost isotropic 3D turbulence. Some of subrotors that do penetrate into the interior of the main rotor, such as s4 in Fig. 7b, are temporarily intensified by vortex stretching in the turbulent 3D flow. As a consequence, the maximum along-ridge vorticity of the subrotors in the 3D simulations is stronger than that in the 2D case.
A number of additional 2D simulations were conducted with varying mountain heights in order to examine the sensitivity of the subrotor structure to the lee wave amplitude. In all of the simulations that contained mountains of sufficient height to produce leeside rotors, vortex rollup occurred along the vortex sheet at the leading edge of the rotor. Larger amplitude lee waves produced more vigorous vortices along the vortex sheet. Circulations within the rotor beneath the primary lee wave were present in all simulations that contained rotors. Overall, these simulations indicate the subrotor structure is quite sensitive to the wave amplitude. For relatively low mountain heights (e.g., 1000 m), which are sufficiently high to induce boundary layer separation and rotors, the vortex sheet did not lift as high off the surface after the separation point and the positive vorticity circulation to the rear of the rotor (e.g., Fig. 6a) was wider and not as deep relative to simulations conducted with higher mountains. Some mountain heights result in lee waves and associated pressure perturbations that are favorable for a greater number of subrotor circulations, with up to five existing simultaneously, while other mountain heights produce rotors that contain as few as two simultaneous subrotors.
4. Three-dimensional topography
A series of simulations were conducted over fully three-dimensional ridges to explore the sensitivity of rotors and subrotors to the ridge aspect ratio and along-ridge terrain variations. A diagnosis of the horizontal and vertical structure of the rotors and subrotors for such flows is presented in this section, along with detailed vorticity diagnostics.
a. Sensitivity of rotors to the length of the ridge
Epifanio and Durran (2001) found that the character of flow past a uniform 3D ridge of finite length is sensitive to the horizontal aspect ratio β (the ratio of cross-stream to streamwise topographic scales). In the nonbreaking regime, they found that the flow in the interior of the ridge approached the 2D limit when β ≥ 10. When breaking waves were present, however, significant differences were found between the 2D and 3D cases for aspect ratios at least as large as 12; these differences were attributed to a marked increase in flow deflection around the ridge associated with the onset of wave breaking. Previous 2D rotor simulations of Doyle and Durran (2002) and Vosper (2004) have demonstrated that a fundamental relationship exists between lee wave amplitude, boundary layer separation, and rotor formation. Furthermore, Doyle and Durran have noted that there is a strong correlation between the strength of the lee-wave-induced pressure gradients and the strength of the reversed flow within the rotor. Recognizing the sensitivity of rotors to mountain wave amplitude, based on the previous studies one would expect a close relationship between the rotor strength and the horizontal aspect ratio of the ridge.
A series of three-dimensional simulations were conducted to explore the sensitivity of the rotor characteristics to β. The topography in these simulations is specified by (12) and (13) except that hm = 1000 m. Since our focus is not on subrotor structure and our parameter space includes very long mountains, a larger horizontal domain is used with three nested grids such that the finest mesh, which contains 187 × 277 grid points, is centered on the mountain and has a horizontal grid increment of 600 m. A total of 90 vertical levels are used with a resolution identical to the simulations discussed above. The simulations are integrated to a time of 3 hours.
The simulation results indicate that the rotor strength and depth are, indeed, sensitive to the horizontal aspect ratio of the ridge. The relationship between the maximum reversed flow within the rotor and the depth of the reverse flow are shown in Fig. 8 as a function of β. The metrics are normalized by the maximum reversed flow and reverse flow depth from three-dimensional simulations for an infinitely long uniform ridge with the same ridge profile. A rapid increase in the strength of the reversed flow and depth of the reversed flow occurs as β is increased from unity (the circular mountain case) to 6. The depth of the reversed flow is almost independent of horizontal aspect ratio and slightly less than that in the infinite ridge limit for β > 6. The strength of the reversed flow also is always weaker than in the infinite ridge limit, and gradually approaches the infinite ridge value as β increases from 9 to 20.
b. Sensitivity of rotors to along-ridge topographic variations
As noted in connection with Fig. 7, stretching of vortex tubes along an axis parallel to the ridge appears to intensify some of the subrotors in the 3D simulation for flow over an infinitely long uniform ridge. Along-ridge topographic variations have the potential to focus such vortex stretching in preferred spatial locations, and this possibility is explored by examining a simulation of the rotors generated by flow over the isolated ridge with a localized peak shown in Fig. 3. Note that our analysis is focused on the nested meshes placed just to the south of the peak.
1) Rotor and subrotor structure and characteristics
Figure 9 shows the time-mean characteristics of the rotor circulation on the finest (60 m) mesh along the vertical cross section A–A′ indicated in Fig. 3. The time-mean fields were computed from model output at a frequency of 20 s over the interval between 3 and 3.75 h. The mean horizontal velocity and potential temperature fields are roughly similar to those obtained for 3D flow over the infinitely long uniform ridge shown in Fig. 5c, although the wave amplitude is higher and the reversed surface flow beneath the rotor (which exceeds 10 m s−1) is stronger. The amplitude is higher because the isolated peak just north of cross section A–A′ raised the maximum height of the topography to 2 km from the 1.5-km value for the uniform ridge. Note that the position of the rotor is also shifted closer to the ridgeline in the isolated-peak simulation, which accounts for most of the difference in the topographic profiles plotted in Figs. 5 and 9.
The field of mean y-component horizontal vorticity, shown in Fig. 9b, is also similar to that for the 3D flow over the infinitely long uniform ridge (Fig. 5d), except that due to the increased wave amplitude, the sheet of positive η that is lifted off the surface by the lee wave rises to a height of 2.2 km at its apex, whereas the vortex sheet rises to only 1.4 km in the lee of the uniform ridge. In addition, the vorticity field within the main rotor shows more finescale structure than that in the uniform ridge case (Fig. 5d) because no averaging is performed along the y coordinate.
The evolution of the y-component vorticity at 1-min intervals along cross section A–A′ is shown in Fig. 10 beginning at 3 h 24 min. As with the time-averaged fields, the vorticity distribution is broadly similar to that obtained in the 3D simulation for flow over the infinitely long ridge. Aside from the difference in the height of the main rotor, already discussed in connection with Fig. 9, the main differences introduced by the isolated peak are 1) a tendency toward greater amplification of the subrotors as they roll up and detach from the main vortex sheet along the upstream side of the lee wave and 2) a tendency for the subrotors to weaken more rapidly downstream of the crest of the lee wave. The intensification of subrotors, such as s5 in Fig. 10, is of particular interest because the resulting vorticities, as large as 0.25 s−1, are the strongest encountered in these simulations. As will be demonstrated shortly, this intensification is produced where the flow is modified by the isolated peak through vortex stretching along an axis parallel to the mountain ridge. Note that the character of the simulated subrotors is qualitatively similar to those captured in lidar observations by Banta et al. (1990), which revealed three small-scale circulations present in the streamwise component of the wind and located along a strong shear zone within a rotor during a downslope wind event. The horizontal scale of the observed circulations was roughly 500–1000 m, which is broadly similar to the scale of the subrotor circulations in Fig. 10.
The velocity perturbations associated with the intense subrotor s5 of Fig. 10b are shown in Fig. 11 by vertical cross sections along A–A′ of the cross-mountain wind component and the vertical velocity. Clearly very strong wind shear is associated with this subrotor. A horizontal wind maximum of 26 m s−1 is located just 1.25 km above a region of −8 m s−1 reversed flow. A couplet in the vertical velocity is also apparent with an updraft of over 22 m s−1 and a downdraft of −7 m s−1 separated horizontally by less than 2 km. The strength of this updraft appears physically reasonable, as it is just slightly larger than the maximum vertical gust recorded by the B-29 during the Jet Stream Project (Fig. 1).
Owing to the presence of the isolated peak, the along-ridge variation in subrotor s5 is highly structured, as illustrated in Figs. 12a–c, which shows fields at z = 1.2 km (the level of s5 at t = 3 h 24 min) on a subdomain of the finest resolution mesh (Δx = 60 m). The cross-mountain wind component (Fig. 12a) shows a rapid transition from strong downslope flow, which is over 30 m s−1 near the left edge of the plot, to weak reversed flow within the main rotor beneath the lee wave. Within the main rotor there are several patches of wind directed back toward the mountain with maximum magnitudes exceeding 8 m s−1. Although the variations in the mountain-parallel wind component (Fig. 12b) within the large-scale rotor are roughly half as strong as those in the cross-mountain component, significant variations are still present over distances of a few kilometers.
The y component of the horizontal vorticity (Fig. 12c) is weakly negative owing to the decrease in u with height above the level of maximum downslope winds on the left edge of the plot, and then transitions to a band of strong positive η where the primary vortex sheet rises through the z = 1.2 km level at the leading edge of the main rotor. The vorticity field inside of the main rotor is rather chaotic with numerous smaller-scale features reminiscent of the horizontal vortex tubes found in simulations by Clark et al. (2000) of clear-air turbulence in wave breaking regions in the upper troposphere. At this time, s5 has a maximum horizontal vorticity of 0.21 s−1, a horizontal dimension of approximately 600 m × 2250 m, and is embedded in the vortex sheet at the leading edge of the rotor (near the center of the plot). One minute later (at 3 h 25 min), this subrotor intensifies to a maximum vorticity of 0.25 s−1, elongates in the y direction to over 3000 m, and rises vertically to near the 1.7-km level due to advection by the lee wave (Fig. 12d). Note the vorticity within s5 has intensified to values considerably greater than that of the ambient vorticity within the vortex sheet below.
2) Vorticity generation
The vorticity evolution of three-dimensional flows is considerably more complex than the analogous two-dimensional case because of the contribution of tilting and stretching processes. In three dimensions, the y-component horizontal vorticity tendency budget, η, for the compressible set of equations can be expressed as
where B is the buoyancy given by B = g(θ − θ)/θ, ξ and ζ, are the x- and z-vorticity components, respectively, and Du and Dw are the subgrid-scale dissipation terms for the u- and w-momentum equations, respectively. The first term on the right-hand side of the equation is the stretching term, the second represents the y-vorticity contribution from tilting of vorticity into the y plane from other components, the third term is the baroclinic contribution, and the fourth term is the vorticity generation due to surface friction and subgrid-scale dissipation processes.
To illustrate the vorticity generation for an intense subrotor, the vorticity budget at t = 3 h 24 min, z = 1.2 km over a subsection of the finest resolution (60 m) mesh centered on s5 is shown in Fig. 13. At this time, subrotor s5 has a maximum η of 0.21 s−1 (Fig. 13a) and, as shown in Figs. 10a,b and 12c,d, it is intensifying while rising at the speed of the background flow along the leading edge of the lee wave. The intensification is also apparent in the field of total Lagrangian vorticity tendency (Fig. 13b), in which the region of maximum intensification is roughly coincident with the center of s5. The values in Fig. 13b are computed by subtracting the advective tendency from 1/Δt times the difference in η between two adjacent time steps. As a check on the closure of the vorticity budget, the sum of all of the individual terms in the Lagrangian vorticity budget is plotted in Fig. 13c; clearly there is close agreement with the values in Fig. 13b. Among the individual contributions to the Lagrangian vorticity budget, the stretching term (Fig. 13d) is by far the largest contributor to the intensification of s5, with tilting (Fig. 13e), friction and dissipation (Fig. 13f), and baroclinic processes (not shown) considerably less important. Farther downstream within the circulation of the main large-scale rotor, the tilting term does become more important as turbulent mixing becomes more three-dimensional and the vorticity in other planes can be easily tilted into the y component.
The significance of vortex stretching for subrotor intensification was explored further by examining the vorticity budget for other strong subrotors in the 3D simulation with a nonuniform finite ridge. The strongest subrotors were identified in the simulation at 20-s intervals over the period between 3 and 3.75 h. Only subrotors above 200 m AGL and below 3000 m were considered to avoid maxima embedded within the surface-based vortex sheet and above the primary larger-scale rotor circulation. The vorticity budget for subrotors of varying intensities is shown in Fig. 14. Each term in the budget was averaged over a 5-point horizontal stencil surrounding the center of the subrotor vorticity maximum. The results indicate that the stretching term is the leading contribution to the instantaneous vorticity tendency in 72% of the subrotors, while the tilting of vorticity from the x and z planes is the dominant generation mechanism in the remaining 28% of the subrotors. The baroclinic term makes only very small contribution to η in the region containing the subrotors. The dissipation term acts as a downgradient sink of positive vorticity. For stronger subrotors, having vorticity maxima of 0.2 s−1 or greater, stretching dominates the budget in a slightly larger fraction of the cases (78%).
As suggested previously, along-ridge terrain variations can introduce regions of enhanced stretching deformation, which subsequently may create favorable regions for subrotor intensification. In the incompressible limit, the stretching term in (14) is proportional to ∂υ/∂y, and positive values of ∂υ/∂y are clearly evident in the vicinity of s5 in Fig. 12b. Here we further explore the influence of stretching associated with the presence of the isolated peak by considering the differences in the stretching term in the η vorticity budgets for the northern and southern halves of the second-finest (180 m) resolution mesh. As shown in Fig. 3, there is much more variability in the terrain along the y coordinate in the northern half of the 180-m resolution mesh than in the southern half. Thus, subrotors present in the northern half of the 180-m subdomain may intensify more rapidly through systematic vortex stretching than those in the southern half domain. To examine this hypothesis further, a temporal mean of the model state variables and the terms in the vorticity tendency equation were accumulated at 20-s time resolution over the period between 3 and 3.75 h. The temporal mean quantities were then averaged in the spanwise direction across each of the two individual half domains. Vertical cross sections of the vortex stretching and potential temperature for the northern and southern half domains are shown in Figs. 15a and 15b, respectively. There is systematic stretching of the vortex sheet at the leading edge of the rotor in both halves of the domain, but the mean vortex stretching tendency for the northern half of the domain is more than a factor of 2 larger than that in the southern half.
The relation between the regions of most intense vortex stretching, the lee waves, and the topography is further illustrated in Fig. 16, which shows the −2 K perturbation potential temperature isoline and contours of a stretching metric superimposed on the terrain contours. The perturbation potential temperature is vertically averaged through the layer between 2 and 3 km and temporally averaged between hours 3 and 3.75 using data from the 540-m resolution mesh. The crest of the first lee wave (the region where the perturbation θ is less than −2 K, labeled 1 in Fig. 16) forms an uninterrupted line along the foot of the lee slope. Farther downstream the second, third, and fourth waves are progressively weaker and more disrupted by the isolated peak; they also show evidence of a diverging “ship bow wave” pattern at the north end of the ridge.
The stretching metric Sη is the product max(η, 0) × max(∂υ/∂y, 0) and is averaged over the same layer and time interval used to compute perturbation θ. This metric captures that portion of the vortex stretching that will intensify ridge-parallel positive vorticity (in the Boussinesq limit). The most intense Sη appears within the first lee wave in a patch just south of the isolated peak. Other regions of intense stretching-induced amplification are found under the crest of the first lee wave downstream and adjacent to the isolated peak. Weaker maxima in Sη are found farther downstream of the peak, north of the peak under the first and second lee waves (where large values of ∂υ/∂y occur near the end of the ridge), and at other scattered locations under the crests of waves 2 through 4.
In general, the stretching-induced amplification Sη is strongest directly under the lee wave crest rather than along the leading edge of the wave, as in Fig. 15. This shift in position is due to the relatively high elevation of the layer used to compute Sη. The 2–3-km layer was selected because it was difficult to accurately compute the vertical average of Sη in layers that intersect the terrain. Note that the contour interval for Sη in Fig. 16 is an order of magnitude smaller than that for the time- and y-averaged stretching in Fig. 15; this is because the stretching is weaker in the 2–3-km layer than lower down along the leading edge of the lee wave and because the computations in Fig. 15 are done at three times finer resolution than those in Fig. 16.
5. Subrotor instability
The instability that develops along the vortex sheet, apparent in Fig. 10, is well correlated with strong horizontal wind shear present along the leading edge of the rotor. Possible candidates for an instability mechanism along the vortex sheet include parallel shear or Kelvin–Helmholtz (KH) type of instabilities. The Miles–Howard condition for KH instability (Miles 1961; Howard 1961) states that a necessary condition for KH instability is that somewhere in the flow the Richardson number
Considering the strong vertical wind shear (Fig. 11) and the weak stability within the large-scale rotor circulation (Fig. 9a), layers of small Ri are certainly possible within the rotor, particularly near the leading edge.
To explore the subrotor instability characteristics further, vertical profiles of the potential temperature, horizontal wind components, y component of horizontal vorticity, and Ri were extracted from the 3D simulation at two different y locations along the lee slope downwind of the flow-separation point. At one location strong subrotors were repeatedly generated, whereas subrotors were weaker or absent at the second location. Beginning at t = 3h 20 min, both profiles were averaged over 10-min time periods, corresponding to three subrotor cycles, to better capture the structure of the mean upstream flow.
The potential temperature profile from the location most favorable for subrotor development (solid line Fig. 17a) has a three-layer structure with a neutral layer below 300 m (AGL), an unstable region in the 300–600-m layer, and strong stratification farther aloft. Strong forward vertical shear of the cross-mountain wind component is present also in the 300–600-m layer (Fig. 17b), which is flanked by two weakly sheared layers. The along-ridge velocity component (Fig. 17b) is weak with little vertical shear. The layer with large ∂u/∂z is associated with a well-defined maximum in the y vorticity of the mean state near z = 550 m with weak vorticity above and below, which is consistent with the mean vertical cross sections of η in Fig. 9. A deep layer of Ri ≲ 0.1 extends from the surface to 600 m with a minimum Ri ∼ −1 at z = 370 m (Fig. 17d). The time-averaged profile from the location with weak subrotor events (dashed line, Fig. 17) exhibits weaker vertical shear and increased stability in the middle layer, and the maximum y vorticity is reduced by a factor of 2 relative to that in the profile most favorable for subrotor development.
A series of high-resolution 2D simulations were conducted to explore the dynamic stability of these profiles further. The model was initialized with the horizontally homogeneous profiles shown in Fig. 17 except that the vortex sheet was raised an additional 600 m to prevent the subrotors from interacting with the lower boundary. Random perturbations of maximum magnitude 0.1 m s−1 were superimposed on the initial cross-mountain wind component. The horizontal and vertical resolutions were both 20 m with 601 grid points in the horizontal and 140 in the vertical. The simulations used a no-slip lower boundary condition and a radiation condition at the upper and lateral boundaries except that a Rayleigh damping region was added at the inflow boundary.
Vertical cross sections of η for these simulations are shown at t = 10 and 20 min in Fig. 18. In all cases, vortex rollup occurs along the primary shear zone in a manner essentially identical to KH billows in a stratified flow (e.g., Klaassen and Peltier 1985; Fritts et al. 1996; Palmer et al. 1996). The only aspect of this flow that differs from classical KH problem is the presence of a statically unstable layer in the region of strongest mean-state η (see Figs. 17a,c). Additional simulations (not shown) performed with uniform static stability across the shear layer consistent with the classical KH instability problem demonstrate that the vortex rollup is almost unchanged (although the maximum vorticities are just slightly stronger when the unstable layer is present).
It is apparent that the vortex sheet rolls up much more slowly when the basic state corresponds to the weak subrotor profile (Figs. 18a,b) than when the basic state matches the profile favorable for strong subrotor development (Figs. 18c,d). In the strong subrotor case, the horizontal distance between subrotor centers is roughly 1.8 km and new vortices break off the main vortex sheet every 2–3 min. Despite the upward slope of the vortex sheet and other details about the background flow in the actual rotor simulations, the frequency of new-vortex rollup, and the horizontal and vertical scales of the billows shown in Fig. 18 are roughly similar to those of the subrotors breaking off the vortex sheet in Figs. 6, 7 and 10. Thus, a KH-like instability does appear to be responsible for subrotor formation. An additional 3D simulation (not shown) using a domain size of 601 × 151 × 140 with identical resolution confirmed that the initial breakdown of the shear layer into quasi-2D vortices occurs on the same time scale as in the 2D simulations. In the 3D simulation, however, those initially 2D vortices continue to break down into smaller-scale turbulence as in previous 3D simulations of secondary instabilities associated with KH billows (e.g., Fritts et al. 1996; Palmer et al. 1996).
We have explored the dynamics and internal structure of mountain-wave-induced rotors through a series of very high resolution 2D and 3D simulations using a large-eddy-simulation version of the nonhydrostatic COAMPS model. Consistent with the previous findings of Doyle and Durran (2002), the main large-scale rotor is produced by boundary layer separation, as a vortex sheet generated by mechanical shear in the boundary layer below strong downslope winds is lifted aloft by the lee wave circulation. In these high-resolution simulations, the ascending vortex sheet becomes unstable and breaks up into in a series of subrotors along the leading edge of the larger-scale rotor circulation.
Marked differences in the internal structure of the main rotor and in the subrotor evolution are apparent between the 2D and 3D simulations. The substructures within the main rotor tend to be larger in horizontal scale and steadier in the 2D simulations. In the 2D case the main rotor consists of a large, relatively strong circulation whose vorticity matches the sign of the vorticity in the boundary layer-shear-generated vortex sheet. As subrotors break off the main vortex sheet, they are typically entrained into this circulation and, as they dissipate, contribute toward the maintenance of the vorticity in the main rotor. In contrast, in the 3D simulations, including those for an infinitely long uniform ridge, the vorticity field in the interior of the main rotor is broken into more chaotic finescale structures and its time-mean value is much weaker. In addition, the vorticity within the subrotors appearing along the boundary between the main rotor and the surrounding lee wave is significantly more intense than in the 2D case. For some of the most intense subrotors in the 3D simulations, the along-ridge component of the vorticity approaches 0.3 s−1, and extreme vertical shears in the horizontal wind and horizontal shears in the vertical wind are both produced by vector velocity differences on the order of 30 m s−1 over a distance of about 1 km.
Vorticity budget calculations indicate that stretching of the along-ridge horizontal vorticity and tilting are important sources for intensification of the vorticity within the strongest 3D subrotors. Along-ridgeline variations in the topography in the form of an isolated peak protruding above an otherwise uniform ridge induced preferred regions of horizontal difluence favorable for vortex tube stretching and subrotor intensification. In particular, the mean stretching tendency for η vorticity in the main vortex sheet just south of the isolated peak is more than a factor of 2 larger than in the region farther south where the ridge is more uniform. More strong subrotors were also found to develop in the region just south of the isolated peak than in an area of equivalent size farther to the south.
The instability of the main vortex sheet at the leading edge of the lee wave appears to be essentially a Kelvin–Helmholtz instability. Additional simulations show that perturbations about a horizontally uniform initial state with shear and static stability profiles characteristic of the structure of the main vortex sheet produce a series of rolled-up vortices similar to the subrotors that develop in the full mountain-wave/rotor simulations.
Preliminary analysis of recently collected observations from Doppler lidars (Weissmann et al. 2006) and the Raman-shifted eye-safe aerosol lidar (REAL) aerosol backscatter lidar (De Wekker et al. 2006) during T-REX suggest the presence of small-scale vortices along the leading edge of the main rotor that are at least qualitatively similar to those appearing in our 3D simulations. We anticipate that further analysis of the T-REX observational dataset will help guide future high-resolution numerical simulations of rotor and subrotor circulations.
Joachim Kuttner, Chris Golaz, and Qingfang Jiang are gratefully acknowledged for helpful discussions. The first author acknowledges support through the Office of Naval Research’s Program Element 0601153N. The support for the second author was provided by the National Science Foundation Grant ATM-0506589. Computational resources were supported in part by a grant of HPC time from the Department of Defense Major Shared Resource Centers, Aberdeen, Maryland, and Wright Patterson Air Force Base, Ohio. COAMPS is a registered trademark of the Naval Research Laboratory.
Corresponding author address: James D. Doyle, Naval Research Laboratory, Marine Meteorology Division, 7 Grace Hopper Avenue, Monterey, CA 93943-5598. Email: email@example.com
This article included in the Terrain-Induced Rotor Experiment (T-Rex) special collection.