Numerical simulations of flow over steep terrain using 11 different nonhydrostatic numerical models are compared and analyzed. A basic benchmark and five other test cases are simulated in a two-dimensional framework using the same initial state, which is based on conditions during Intensive Observation Period (IOP) 6 of the Terrain-Induced Rotor Experiment (T-REX), in which intense mountain-wave activity was observed. All of the models use an identical horizontal resolution of 1 km and the same vertical resolution. The six simulated test cases use various terrain heights: a 100-m bell-shaped hill, a 1000-m idealized ridge that is steeper on the lee slope, a 2500-m ridge with the same terrain shape, and a cross-Sierra terrain profile. The models are tested with both free-slip and no-slip lower boundary conditions.
The results indicate a surprisingly diverse spectrum of simulated mountain-wave characteristics including lee waves, hydraulic-like jump features, and gravity wave breaking. The vertical velocity standard deviation is twice as large in the free-slip experiments relative to the no-slip simulations. Nevertheless, the no-slip simulations also exhibit considerable variations in the wave characteristics. The results imply relatively low predictability of key characteristics of topographically forced flows such as the strength of downslope winds and stratospheric wave breaking. The vertical flux of horizontal momentum, which is a domain-integrated quantity, exhibits considerable spread among the models, particularly for the experiments with the 2500-m ridge and Sierra terrain. The differences among the various model simulations, all initialized with identical initial states, suggest that model dynamical cores may be an important component of diversity for the design of mesoscale ensemble systems for topographically forced flows. The intermodel differences are significantly larger than sensitivity experiments within a single modeling system.
The fundamental linear theory for the generation of inviscid mountain waves forced by stratified airflow over two-dimensional obstacles has been established for several decades (e.g., Queney et al. 1960; Smith 1979, 1989). Vertically propagating mountain waves often amplify in the stratosphere due to the decrease of atmospheric density with altitude and nonlinear processes, which may lead to overturning and turbulent breakdown (e.g., Lindzen 1988; Bacmeister and Schoeberl 1989). Mountain waves can have an important impact on the atmosphere because of their role in downslope windstorms (Klemp and Lilly 1975); clear-air turbulence (Clark et al. 2000); vertical mixing of water vapor, aerosols, and chemical constituents in the stratosphere (Dörnbrack and Dürbeck 1998); potential vorticity generation (Schär and Durran 1997); and orographic drag influence on the general circulation (Bretherton 1969; Ólafsson and Bougeault 1996).
Although numerical models have been able to successfully simulate gravity wave characteristics (e.g., Smith and Grønås 1993; Schär and Durran 1997) including breaking (Clark and Peltier 1977; Bacmeister and Schoeberl 1989), and other related topographically forced flows (e.g., Sharp and Mass 2004; Colle and Mass 1996; Doyle et al. 2005; Colle et al. 2008), the degree to which mountain waves are predictable still remains an outstanding question. Lorenz (1969) suggested that perturbation growth rates increase as the horizontal scale of the phenomenon decreases, effectively limiting the intrinsic predictability. The notion that mesoscale phenomena forced by the lower boundary attain enhanced predictability has been hypothesized (e.g., Anthes et al. 1985). However, this perspective is likely overly optimistic because of lateral boundary conditions, numerical dissipation, and adjustment issues (Errico and Baumhefner 1987; Vukicevic and Errico 1990), as well as nonlinearities introduced by the underlying terrain. For example, mesoscale predictions of landfalling fronts were found to be very sensitive to small changes in incident flow, as deduced through simulations made with small modifications to the topography orientation by Nuss and Miller (2001). Two-dimensional idealized adjoint (Doyle et al. 2007) and ensemble (Doyle and Reynolds 2008) model results indicate large sensitivity to the initial state as the mountain height increases, forcing wave breaking, where perturbation growth becomes extremely rapid leading to large ensemble spread. Reinecke and Durran (2009a) used three-dimensional ensemble simulations to assess the sensitivity of downslope winds to the initial conditions. They examined two downslope windstorm cases and found the results to be quite sensitive to the initial state, especially for an event that featured low-level wave breaking.
An intercomparison of model simulations of the 11 January 1972 Boulder downslope windstorm is discussed in Doyle et al. (2000). Upper-level wave breaking was predicted by all of the models in comparable locations. However, there were a number of significant differences among the simulations including the details of the wave breaking and the lower-tropospheric wave structure, which was characterized by a hydraulic jump in most of the models and large-amplitude waves in several of the simulations. In the Doyle et al. (2000) study, only a single model test case was carried out. A number of outstanding questions remain including sensitivity of mountain-wave predictions to the model formulation.
During the Terrain-Induced Rotor Experiment (T-REX; Grubišić et al. 2008), high-resolution forecasts were routinely conducted to assist in mission planning using a number of different three-dimensional nonhydrostatic numerical models such as the Coupled Ocean–Atmosphere Mesoscale Prediction System (COAMPS1; Hodur 1997), two dynamical cores of the Weather Research and Forecasting model (WRF), namely the Advanced Research WRF (WRF-ARW; Skamarock et al. 2008) and the Nonhydrostatic Mesoscale Model (NMM) (Black et al. 2005; Janjic et al. 2010), as well as the fifth-generation Pennsylvania State University–National Center for Atmospheric Research (PSU–NCAR) Mesoscale Model (MM5; Dudhia 1993). The forecasts of topographically forced flows, such as windstorms and mountain waves, often differed significantly among the various three-dimensional models. Although one should not necessarily expect forecasts from multiple models to be in close agreement (e.g., Dörnbrack et al. 2005), these differences among the models motivated a more fundamental intercomparison of models in a more controlled and idealized setting. In this study, we compare the results of two-dimensional numerical simulations using an observed initial state from T-REX in order to further investigate the sensitivity of downslope windstorm and mountain-wave predictions to model characteristics through a series of model test cases. Strong downslope windstorms and complex regions of upper-level wave breaking develop in some of the test cases. The evolution of such phenomena presents a particularly difficult challenge from a mesoscale predictability perspective. The overall goal of this study is to explore the predictability and characteristics of model solutions for conditions under which wave breaking and a downslope windstorm were both observed and anticipated from theoretical considerations. The experimental design and model descriptions are presented in section 2. The simulation results are described in section 3, followed by the discussion and conclusions in section 4.
2. Experimental design and description of numerical models
In these tests, the configuration of the individual models is standardized as much as possible to facilitate straightforward comparisons and analysis. The models are applied in a two-dimensional mode with a horizontal grid increment of 1 km and a vertical grid spacing that is stretched from 50 m near the surface to a constant of 200 m, which extends from an altitude of 800 m to at least 26 km (134 or more vertical levels). The fine vertical grid increment in the lowest 800 m of the model is used to better resolve the near-surface processes, particularly for the simulations performed with surface friction.
The attributes of the 11 numerical models that are used in this intercomparison study are summarized in Table 1. The numerical model suite includes: the Advanced Regional Prediction System (ARPS; Xue et al. 2000, 2001), All-Scale Atmospheric Model (ASAM; Hinneburg and Knoth 2005), Boundary Layer Above Stationary, Inhomogeneous Uneven Surfaces (BLASIUS; Wood and Mason 1993), two different versions of COAMPS (versions 3 and 4; Hodur 1997), the Durran and Klemp (1983) model (DK), the Bryan cloud model (CM1; Bryan and Fritsch 2002), the Eulerian/semi-Lagrangian model (EULAG; Prusa et al. 2008), the Regional Atmospheric Modeling System (RAMS; Pielke et al. 1992), the Met Office Unified Model (UM; Davies et al. 2005), and the WRF-ARW model (Skamarock et al. 2005). These models are quite diverse with regard to their numerical characteristics. For example, the suite of models includes semi-Lagrangian (UM) and Eulerian models, models with higher-order time differencing (CM1, WRF), anelastic (BLASIUS, EULAG) and compressible (e.g., DK, CM1) dynamical cores, and vertical coordinates with one model that uses a cut-cell technique (ASAM) and the others that follow the terrain. The models have a diverse set of turbulence parameterizations, as well as treatment of the lateral and upper-boundary conditions. Simulations are conducted with zero surface stress and a bulk stress parameterization, which we refer to as free slip and no slip, respectively. No rotation is used for the simulations with a free-slip lower boundary. A Coriolis parameter consistent with a 37°N latitude is used for the experiments with surface friction. The Coriolis force is applied to the perturbation velocity rather than the mean fields in the no-slip simulations, hence the basic-state flow is assumed to be in geostrophic balance. The surface heat flux is zero in all models, implying that the ground temperature is in balance with the surface air temperature. No moist processes are included in any of the simulations. Vertical mixing in the free atmosphere is the only physical parameterization included within the models for the free-slip simulations. For the no-slip simulations, the roughness length is specified at 10 cm for most of the models, with the exception of EULAG and CM1, which use bulk aerodynamic drag coefficients of Cd = 10−3 and 0.008, respectively.
A total of 6 test cases were conducted using the 11 different numerical models. The setup parameters for the test cases are contained in Table 2. The first test case consists of a baseline experiment for hydrostatic mountain waves that uses a constant isothermal temperature profile of 250 K and a wind speed of 20 m s−1. The mountain profile in the idealized topography cases is a two-sided Witch of Agnesi profile:
with a = au for x < 0 and a = ad otherwise. In the baseline test case, the mountain is symmetric with a = 10 km. The terrain shape we consider for the Ex1000_fs, EX2500_fs, and EX2500_ns cases is asymmetric with the lee slopes steeper than the upwind slopes, which is a characteristic of a number of mountain ranges including the Colorado Front Range and the Sierra Nevada. To represent this asymmetry in these tests, the upwind half-width is set to au = 40 km and the downwind half-width to ad = 5 km (following Hertenstein and Kuettner 2005). In the baseline case, the mountain height is hm = 100 m, which generates small-amplitude mountain waves in the linear hydrostatic regime for the reference sounding conditions considered. The Ex1000_fs case makes use of a 1000-m maximum mountain height, and hm = 2500 m is used in the Ex2500_fs and Ex2500_ns experiments. The terrain used in ExSierra_fs and ExSierra_ns is based on a section across the Sierra Nevada and oriented along the T-REX “B” flight track (e.g., Grubišić et al. 2008; Smith et al. 2008). A 1-km resolution terrain digital elevation model [the National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center Global Land One-kilometer Base Elevation (GLOBE)] is used to specify the terrain with a maximum height of ~3500 m. The terrain upstream of the Sierra Nevada and downstream of the second mountain range, the Inyos, is eliminated to isolate the local Sierra response. The terrain used in the ExSierra_fs and ExSierra_ns simulations is filtered to remove 2Δx variations.
The initial conditions for all experiments are horizontally homogeneous and with the exception of the baseline are based on the 2100 UTC 25 March 2006 Mobile GPS Advanced Upper-Air Sounding System (MGAUS) sounding upstream of the Sierra Nevada during T-REX IOP 6 (Fig. 1), which was one of the stronger mountain-wave cases observed during T-REX (Grubišić et al. 2008; Smith et al. 2008; Doyle et al. 2011). The wind component parallel to the T-REX “B” flight track (245°) is used for the initial state. The MGAUS sounding is extended in the vertical through the use of the National Weather Service (NWS) Oakland radiosonde profile above 14 km. The wind profile was smoothed in the low levels to remove several sharp shear layers and a minimum wind speed of 5 m s−1 is used near the surface. Near the top of the sounding above 18 km, the wind speed decreases linearly with height to 0 at 23 km to minimize problems with gravity wave reflection from the top boundary. The sounding contains multiple shear layers (Fig. 1), strong low-level cross-mountain winds (>20 m s−1) for the simulations that have terrain heights at 2500 m or above, and increased stability above the mountain crest between 1 and 4 km (Figs. 1 and 2), with the latter two characteristics favorable for downslope windstorms (e.g., Brinkman 1974; Durran 1986, 1990). As discussed above, a critical level is present at ~23 km where the cross-mountain wind speed is small. The Scorer parameter, , where N is the Brunt–Väisälä frequency; U is the cross-mountain wind speed; and Uzz is the second derivative of U with respect to height z, which exhibits a rapid decrease with height between 2–5 km (Fig. 1) consistent with the possibility of trapping of low-level wave energy (Scorer 1949). The average Scorer parameter in the upper troposphere (i.e., 5–10 km) is approximately 0.000 35 m−1, implying that only hydrostatic waves with wavelengths longer than the cutoff wavelength, , can penetrate this layer.
3. Simulation results
a. Baseline experiment
The purpose of the baseline experiment is to provide an assessment of the model dynamics through the simulation of a common test case in which an analytic solution can be found. In this case, we selected a linear hydrostatic gravity wave test case, which nearly all models have used previously as a benchmark during their various developmental stages. The vertical velocity at the 4-h simulation time for the baseline experiment using hm = 100 m is shown in Fig. 3 for all 11 numerical models. In general, all of the simulations feature vertical propagating hydrostatic gravity waves similar to the linear analytic solution. The vertical wavelength for all of the simulations is approximately 6.4 km, in agreement with expectations from linear theory, λz = 2πU/N, where N is the Brunt–Väisälä frequency, which for an isothermal, dry atmosphere can be expressed as
The reference state of and U = 20 m s−1 yields N = 0.0196 s−1 and λz = 6.4 km.
Overall the simulation results are in general agreement amongst the models; however, there are several notable exceptions. The BLASIUS model appears to excessively amplify the vertical propagating wave with altitude, a consequence of the model equation set formulation (see below). The RAMS model simulation shows some disagreement with the analytic solution near the top and bottom portion of the sloped ascent and descent regions, which may be arising from wave reflections from the top boundary. The other model simulations in general are in close agreement with the analytic solution.
The relatively rapid growth of wave amplitude with height in the BLASIUS simulations can be shown to be due to the formulation of the anelastic equation set used in this model. The pressure gradient term in the anelastic form of the BLASIUS momentum equations is written as [see Eq. (A1) in the appendix], where p′ is the perturbation pressure (to a reference, hydrostatically balanced profile) and is a reference density, which is a function of height only. The consequence of the different forms of the pressure-gradient term in the momentum equation becomes clear when linearizing the equations and deriving the vertical-structure equation for small-amplitude mountain waves. The BLASIUS equations result in a fundamentally different form for the vertical structure equation, solutions to which grow relatively rapidly with height. The wave amplitude has a dependence on , where ρ0 is a constant reference density. In contrast, solutions to the vertical structure equation derived from either the Lipps and Helmer equations, or Eqs. (A1)–(A3) when the pressure gradient is expressed as , have a dependence on and hence grow less rapidly with height. Further details are provided in the appendix. Linear solutions (not shown) for the constant background flow baseline case derived with the BLASIUS pressure-gradient term agree closely with the BLASIUS results in Fig. 3, whereas solutions based on either or the Lipps and Helmer term agree closely with the results produced by all other models.
It is clear that the cause of the larger wave amplitude at upper levels in the BLASIUS simulation is a direct consequence of the form of the pressure-gradient force used when the model is run in anelastic mode. Note, however, that the BLASIUS model is rarely run in this configuration. Previous mountain-wave studies based on BLASIUS simulations have all used the Boussinesq formulation of the model, in which is assumed independent of height (e.g., Wells et al. 2008; Vosper 2004). The issue does not apply to these studies. The fact that the EULAG results for the baseline case are closer to those of the compressible models and the analytic solution, is consistent with this analysis (see Fig. 3).
b. Ex1000_fs simulations
The simulated potential temperature and vertical velocity at 4 h are shown in Fig. 4 for the Ex1000_fs experiment, which use hm = 1000 m and a free-slip lower boundary. The models in general exhibit a series of trapped waves in the lower troposphere with a wavelength of approximately 10 km. The trapped waves result from wave energy being partially reflected or ducted as a result of the aforementioned decrease in the Scorer parameter with altitude. The number of large-amplitude wave crests present in each model simulation varies considerably. Most models have at least 3 or 4 larger waves, and some models have as many as 5 or 6 waves, particularly for those models using higher-order numerical methods (e.g., WRF-ARW and CM1). The RAMS and ASAM models show considerable damping, evident from a comparison of maximum surface wind speeds in Table 3 and relative to the mean of all of the models (Fig. 4, lower right). Because of their relatively short wavelengths, the trapped waves are only just adequately resolved using the 1-km horizontal grid spacing and thus may be highly sensitive to the numerical methods applied in the models.
At upper levels, vertically propagating waves that leak through the wave duct are apparent in all of the models. The tilt of the wave phase lines with altitude increases in the stratosphere in all models associated with the increase in static stability. The amplitudes of the waves vary considerably once again, with the weakest wave activity found in ASAM and the most extensive regions in the BLASIUS simulation. Most of the models exhibit localized regions of wave breaking in the stratosphere just above the mountain or downstream of it; however, the amplitude of the breaking varies considerably among the models.
c. Ex2500_fs simulations
The potential temperature and horizontal wind speed perturbation for the 4-h simulation time are shown in Fig. 5 for the Ex2500_fs experiment, which utilizes hm = 2500 m and a free-slip lower boundary. With the larger mountain height, the model simulations have considerably larger-amplitude waves relative to the Ex1000_fs simulation. A strong response is apparent in nearly all of the models. Eight of the models, ARPS, BLASIUS, COAMPSv3, COAMPSv4, CM1, DK, EULAG, and WRF-ARW, all indicate the presence of severe winds in the lee of the terrain with maximum near-surface winds in excess of 60 m s−1. These models indicate the presence of low-level breaking, a weak or reversed wind speed minimum in the midtroposphere above or in the lee of the mountain crest associated with the breaking, trapped or secondary waves in the lower stratosphere, and vertically propagating waves and overturning in the stratosphere. It is interesting to note that COAMPSv4 has one of the stronger windstorms in this group of models, and COAMPSv3 has the weakest windstorm in this group of 8 models. The BLASIUS model has a bore or hydraulic-like jump feature that propagates significantly downstream of the mountain, in contrast to the COAMPSv3 windstorm that is confined to the lee slope.
A group of three models have a quite weak response relative to the other simulations (see Table 3) and the mean (Fig. 5, lower right). These models, ASAM, RAMS, and UM, indicate some enhancement of the winds along the lee slope, but differ substantially in character relative to the other simulations. Both the UM and RAMS simulations have a wind speed maximum near the mountain crest. Interestingly, the mountain-wave response in the ASAM simulation is the weakest in amplitude and contains a small-amplitude lee-wave train.
The characteristics of the stratospheric wave breaking vary considerably among the simulations. All of the models have signatures of stratospheric vertically propagating waves. The models with the strongest lower-tropospheric waves tend to produce the largest wave amplitudes in the stratosphere. Nearly all of the models are in agreement with respect to the vertical wavelength, tilt of the wave phase lines, and preferred regions of wave amplification and overturning in the stratosphere. All of the models with the exception of the ASAM and UM simulate secondary wave generation downstream from the wave-breaking region in the lower stratosphere.
The variations among the simulations are further illustrated in the vertical velocity and potential temperature fields at the 4-h time for Ex2500_fs shown in Fig. 6. A number of models that simulate a strong mountain-wave response (e.g., COAMPSv3, COAMPSv4, CM1, DK, EULAG, and WRF-ARW) have global vertical velocity maxima well downstream of the terrain crest, while the ARPS, RAMS, and UM simulations tend to have maxima near or immediately above the terrain crest, more similar to a hydrostatic wave.
d. Ex2500_ns simulations
The u component of the horizontal wind speed perturbation and potential temperature at the 4-h simulation time are shown in Fig. 7 for the Ex2500_ns experiment, which corresponds to hm = 2500 m, a no-slip lower boundary condition and rotation. In general, the models are more in agreement with each other (Table 3) than in the corresponding free-slip simulation (Ex2500_fs). A comparison of the two experiments underscores the profound impact of surface friction and rotation on the mountain-wave response. In general, trapped lee waves are present in all of the simulations, with the number of simulated crests varying among the models, although to a lesser extent than in Ex1000_fs. The exceptions are ARPS and BLASIUS, both of which simulate a single larger-amplitude wave crest, and exhibit a downstream wave train in the lower troposphere that is heavily damped, in contrast to the robust trapped waves apparent in the other simulations. The WRF-ARW exhibits a jetlike structure in the perturbation wind field that extends from above the mountain to the upstream lateral boundary at an approximate altitude of 5 km. The jet results from an in situ increase in momentum at the upstream boundary that subsequently self-advects toward the ridge. The jet has a minimal impact on the solution by 4 h.
The stratospheric wave characteristics are broadly similar in Ex2500_ns; however, there are several distinct differences among the simulations. Several of the models (e.g., BLASIUS, COAMPSv3, COAMPSv4, DK, UM, and WRF-ARW) all contain a large-amplitude and short horizontal wavelength feature apparent in the isentropes in the lower stratosphere positioned above the lee-wave crests. All of the models contain steepening of the isentropes in the stratosphere associated with breaking of the vertically propagating waves; however, the spatial extent and depth of the breaking layer vary substantially. The common stratospheric wave characteristics are apparent in the mean.
e. ExSierra_fs simulations
The ExSierra_fs experiment makes use of a terrain transect across the Sierra and a free-slip lower boundary condition. The simulated potential temperature and horizontal wind speed perturbation for the 4-h time are shown in Fig. 8. Broad similarities are apparent between all of the simulations; for example, they contain significant wave activity in the lee of the Inyo Range, as well as large-amplitude waves in the stratosphere. However, substantial variations are also apparent among the 11 model simulations. The COAMPSv3, COAMPSv4, DK, and EULAG simulations all develop a strong windstorm along the lee slope of the Sierra and confined to the Owens Valley. The BLASIUS simulation exhibits shooting flow along the Sierra lee slope that extends over the Inyo Range and farther downstream, evolving abruptly into a hydraulic-like jump feature. It is interesting to note that the CM1 and WRF-ARW simulations were among the strongest downslope windstorms in Ex2500_fs, while in ExSierra_fs the CM1 and WRF-ARW-simulated wave responses are qualitatively among the weakest using the same basic model setup and diffusion coefficients. Several simulations contain relatively small-amplitude trapped lee waves downstream of the Inyos (e.g., ASAM, UM, and WRF-ARW) in contrast to the models that exhibit strong low-level wave breaking in the lee of the Sierra (DK) and Inyos (BLASIUS). The ASAM model simulation for Ex2500_fs is one of the weakest wave responses, while the ASAM simulation for ExSierra_fs is qualitatively similar to other models such as ARPS, COAMPSv3, EULAG, RAMS, UM, and WRF-ARW (summarized in Table 3).
The stratospheric wave activity varies substantially in the models. Several models, including COAMPSv4, DK, and WRF-ARW, simulate vigorous wave breaking in the stratosphere, in contrast to other models that feature smaller-amplitude waves in the lower troposphere. Once again the basic properties of the vertically propagating stratospheric waves are broadly similar (cf. the mean), with large differences in the finescale aspects, comparable to the Ex2500 experiments.
f. ExSierra_ns simulations
The ExSierra_ns experiment uses an identical Sierra terrain transect as the previous experiment, a no-slip lower boundary condition and rotation. The simulated potential temperature and u component of the horizontal wind speed perturbation for the 4-h time are shown in Fig. 9. Once again the introduction of surface friction and rotation results in a significant reduction in the gravity wave response in all the model simulations, similar to the pair of Ex2500 experiments. The models are in much closer agreement with each other (Table 3 and mean in Fig. 9) relative to the corresponding free-slip experiment (Fig. 8). The ARPS, ASAM, COAMPSv4, CM1, DK, and WRF-ARW simulations exhibit well-defined trapped waves in the lee of the Inyo ridge. In contrast, other models such as BLASIUS, COAMPSV3, EULAG, RAMS, and UM exhibit much less wave activity downstream of the Inyos.
The stratospheric wave characteristics such as wavelength and amplitude (as diagnosed from the potential temperature perturbations) are overall quite similar, although the details of the wave fields differ. All of the models exhibit a large potential temperature perturbation in the 11–13-km layer downstream from the Sierra crest above Owens Valley or the Inyo Range. Similar regions of wave breaking above 15 km are present in all of the models as well.
4. Discussion and conclusions
In addition to analyzing the instantaneous flow fields in the simulations, it is also illustrative to consider their time evolution. The simulated potential temperature and horizontal wind speed for the 1–8-h times are shown in Fig. 10 for the DK model for the Ex2500_fs experiment. Near the early portion of the simulation, the lee-side windstorm is substantially weaker (less than 50 m s−1 at the 2-h time) and the region of wave breaking first develops in the stratosphere. The 3-h simulation time shows wave breaking within the troposphere with a deep region of wave overturning and a stronger lee-side wind response. A hydraulic-like transition occurs when the fast shooting flow rapidly decelerates downstream, particularly evident beginning at 4 h. This feature, similar to a hydraulic jump, continues to propagate downstream with time and is located approximately 100 km downstream from the mountain crest by 7 h. At no point does the simulation revert back into a weaker windstorm state with an absence of tropospheric breaking (e.g., 2-h time) as in the RAMS and UM simulations (Fig. 8). Although the temporal evolution likely does vary between the various model simulations, it appears that these temporal variations are not rapid or large enough to explain the large variability among the model simulations. Some of the differences in the wave amplitudes may arise from the numerical accuracy used in the dynamical cores (Reinecke and Durran 2009b), although the wave strength and accuracy of the numerical schemes are not obviously correlated for these experiments, perhaps due to the contribution of other factors such as diffusion.
The vertical velocity standard deviation of the model simulations for each experiment is shown in Fig. 11. The smallest standard deviation is the baseline experiment, with the greatest variability within the vertically propagating gravity wave above 10 km and positioned over the mountain. The standard deviation for the baseline experiment is approximately two orders of magnitude smaller than that of any other experiment, which reinforces that the models are generally effective at replicating the basic characteristics of the hydrostatic wave case. In Ex1000_fs, the standard deviation is a maximum in the lower troposphere in the lee of the 1000-m mountain because of the phase and amplitude differences in the simulated trapped wave train.
The standard deviation is the largest in the Ex2500_fs and ExSierra_fs simulations, which contain the largest amplitude waves. In Ex2500_fs, the vertical velocity variance is a maximum in two regions: in the midtroposphere in the breaking region and in the stratosphere downstream of the breaking in the zone that contains the secondary short wavelength gravity waves. The ExSierra_fs exhibits two variance maxima: one near and just downstream of the Sierra in the lower troposphere and a second considerably farther downstream beyond the Inyo Range, associated with the trapped waves that generally differ in character (e.g., wavelength and amplitude) among the models. Considering that the wave amplitude increases with the mountain height, we normalize the standard deviations in Fig. 11 by the corresponding mountain height (note: according to linear theory, the wave amplitude is proportional to the mountain height). The resultant normalized deviations (not shown) are still significantly larger in the higher-terrain simulations, indicating that the spread increases with the mountain height much faster than the wave amplitude, presumably due to nonlinear processes.
The two no-slip simulations with large terrain, Ex2500_ns and ExSierra_ns, exhibit less than half of the vertical velocity standard deviation from that of the free-slip companion simulations. The maximum vertical velocity standard deviation is in the midtroposphere and lower stratosphere just above the lee slopes in Ex2500_ns. In contrast, the maximum standard deviation in ExSierra_ns is primarily confined to the upper troposphere and lower stratosphere, where differences in the characteristics of the wave amplification and breaking among the model simulations are the dominant signatures. It is also noteworthy that the vertical velocity standard deviation maximum seems to shift to higher altitudes as a result of surface friction in ExSierra_ns, as well as Ex2500_ns, underscoring the fundamental change in the flow imposed by the no-slip lower boundary.
The representation of the vertical flux of horizontal momentum due to orographic effects is crucial for the skillful prediction of the large-scale general circulation (e.g., Palmer et al. 1986; Lott 1995, Kim et al. 2003). The momentum flux is a metric of the overall gravity wave activity, and is computed as follows:
where the primes are deviations from the two-dimensional domain average; and u and w are the horizontal and vertical wind components, respectively. Vertical profiles for each model simulation are shown in Fig. 12, corresponding to the Ex1000_fs, Ex2500_fs, Ex2500_ns, ExSierra_fs, and ExSierra_ns simulations. The profiles are generally negative with the largest magnitudes near the surface, in line with theoretical considerations. The larger near-surface negative values for the larger mountain heights are consistent with the expectation of a net drag that the topography imparts on the westerly flow.
Overall, the model results for Ex1000_fs (Fig. 12a) show less spread than in other cases, particularly if BLASIUS and COAMPSv3 are excluded. The other cases exhibit considerably larger variations among the momentum flux diagnostics. The largest magnitude of the momentum flux occurs for the high mountain free-slip simulations (Figs. 12b,d). The negative momentum fluxes are considerably larger for the free-slip simulations (Figs. 12b,d) than for the corresponding no-slip simulations (Figs. 12c,e), in qualitative agreement with previous modeling studies (e.g., Ólafsson and Bougeault 1997). The spread among the models is larger for ExSierra_ns relative to Ex2500_ns. The BLASIUS model appears to exhibit a momentum flux minimum near the tropopause and lower stratosphere for both Ex2500_ns and ExSierra_ns. The large spread in the momentum flux among the models initialized from identical initial states supports the application of stochastic parameterization approaches (Palmer 2001) to gravity wave drag representation in large-scale weather and climate models to better represent model uncertainty (Doyle and Reynolds 2008).
To better assess the time-dependent nature of the simulations, we have analyzed one set of simulations (Ex2500_fs) in greater detail based on a subset of 6 models (see legend in Fig. 13) using a 10-min time resolution. The higher time resolution output is only available for the six model subset. Figure 13 shows the characteristics of the subset mean and standard deviation of the vertical momentum flux for both the subset and entire model set. The vertical momentum flux averaged over the 4–5-h time period for each of the 6 models (Fig. 13a) exhibits considerably less spread than is apparent at the 4-h time for the full model set (Fig. 12b). This result suggests that the temporal variability may be an important component of the overall model spread evaluated at a particular time. The momentum flux standard deviation for the 6-model subset, based on the average of each model over the 4–5-h time (green curve in Fig. 13b), is actually quite small, which once again underscores the significance of the transience within the subset models. The intermodel standard deviation of the momentum flux at the 4-h time is similar to the standard deviation due to transience over the 4–5-h time period for the subset (red and blue curves Fig. 13b). However, it should be noted that several of the models not included in the subset have a relatively small-amplitude gravity wave response for Ex2500_fs (e.g., ASAM, RAMS, and UM in Fig. 5). The temporal variability for these models is the same size or more likely smaller than the other larger-amplitude models contained in the subset. Thus, given the diversity of the gravity wave response between models within the full model set, transience likely explains only a portion of the total variance. In support of this notion, the standard deviation of the vertical momentum flux for the 4–5-h time period, averaged over the subset (red curve Fig. 13b), is substantially smaller than the standard deviation computed over the entire model set at the 4-h time (black curve Fig. 13b). This suggests that although the transience does contribute to the variance, there is considerable variance among the model simulations due to the properties of the models.
To further elucidate the sensitivity of the model simulations to the turbulence parameterization and the model dynamics, additional experiments have been performed using EULAG for the same setup as Ex2500_fs. The horizontal perturbation wind components for these tests are shown in Fig. 14 at the 4-h time. The control simulation is shown in Fig. 14a. Simulations were conducted with a Smagorinsky closure (Fig. 14b) and no explicit turbulence closure [implicit large eddy simulation (LES); Fig. 14c]. The results show that the simulations in this experiment are not particularly sensitive to the turbulence parameterization. There are some details in the flow that are different, such as the downstream penetration of the windstorm (cf. Figs. 14a–c) and the structure within the wave-breaking layers. However, the overall character of the key simulation features (e.g., windstorm and wave-breaking strength, depth, upper-level breaking location, downstream response) are quite similar. A simulation conducted with WRF-ARW using a Smagorinsky closure confirmed the relative insensitivity of the simulated gravity waves to the turbulence parameterization for Ex2500_fs (not shown).
EULAG simulations conducted with a deeper upper-absorbing layer (Fig. 14d) that starts at z = 22 km, compared to the control simulation z = 25 km show that the model simulation is also relatively insensitive to the sponge layer formulation. Additional EULAG experiments using a different pressure solver and two different methods to specify the mean state potential temperature formulation yield nearly identical results as the control simulation. These results suggest that the intermodel differences are significantly larger than simulations conducted with various physics and dynamics options exercised within a single modeling system.
The simulations carried out in this study provide a perspective on the predictability of mountain waves, wave breaking, and downslope windstorms using a suite of high-resolution nonhydrostatic models applied in a two-dimensional configuration. These test cases illustrate that as the mountain height is increased beyond the 100-m baseline and 1000-m Ex1000_fs experiments, the flow response increases considerably in complexity and the model diversity becomes more significant. For example, the higher mountain experiments exhibit a strong lee-side windstorm that is present with a deep internal hydraulic-like jump downstream of the terrain. Wave breaking forms above the mountain lee slopes in the midtroposphere and in numerous layers in the stratosphere that spread perturbation energy in the upstream and downstream directions. Secondary gravity waves are generated by the wave breaking and propagate along the tropopause interface. The source of the trapping of the wave energy aloft is the vertical variation of static stability induced primarily by the stratospheric wave breaking, which is a fundamentally nonlinear phenomenon. The spread among the model simulations is quite large for these high mountain experiments, particularly for the free-slip lower boundary condition cases which exhibit vertical velocity standard deviations that are approximately twice as large as the companion no-slip experiments.
It is difficult to identify the key components of models that explain the differences among the simulations, particularly for the high-terrain experiments (e.g., Ex2500 and ExSierra). The model formulations (Table 1) differ substantially, especially related to the model numerical techniques used for time differencing and advection, as well as the vertical diffusion and boundary layer mixing. The fact that the free-slip experiments exhibit the largest differences among the simulations, suggests that the numerical techniques used in the diverse set of models play a dominant role in determining the variability, rather than the model physical parameterizations. The possibility does exist that some of the results may be influenced by specifics of the model setup (e.g., slight differences in surface drag formulations or the lateral boundary conditions); however, the model configurations were standardized as much as possible to minimize inconsistencies.
The simplifying assumptions of two-dimensionality considered here may limit the generality of the results. Future studies of mountain waves and topographic circulations should be conducted in three dimensions with realistic boundary layers to better understand the predictability implications.
Quantifying mesoscale predictability for terrain-forced flows is a formidable challenge. The results of this study build on previous two-dimensional topographical flow ensemble simulations (e.g., Doyle and Reynolds 2008) that suggest that the predictability of wave breaking, downslope windstorms, and mountain-wave-induced turbulence is limited, particularly when wave steepening or breaking and the development of shooting flow down lee slopes are essential ingredients. In this study, the 11 numerical models provide additional insight into the contribution of model error. For relatively low mountain heights in the linear regime, the spread among the simulations is small. However, the spread increases substantially for gravity waves forced by high mountains in the nonlinear regime, underscoring the uncertainties introduced by different formulations of models. Given the sensitivities to the initial state for terrain-forced flows (Doyle and Reynolds 2008; Reinecke and Durran 2009a), as well as the contributions from the model formulation, prediction of the location and timing of wave breaking and associated turbulence, which are of importance for aviation, may only be possible through probabilistic approaches.
We acknowledge the contribution of the T-REX scientists, forecasters, staff, NCAR staff, and flight crews. The first three authors acknowledge support through the Office of Naval Research’s Program Element 0601153N. J. Schmidli was supported by the Swiss National Science Foundation Grant PA002-111427. ARPS is developed by the Center for Analysis and Prediction of Storms, University of Oklahoma. Brian Jamison provided invaluable assistance with the WRF-ARW runs. In the early stages of this work, Prof. Grubišić was supported by an NSF Grant ATM-0524891 to the Desert Research Institute. Ivana Stiperski acknowledges support of the Croatian Ministry of Science through Grant 004-1193086-3036 to the Croatian Meteorological and Hydrological Service. S. Zhong acknowledges the support of the National Science Foundation Grant ATM-0646299. Drs. Kirshbaum and Zhong’s simulations were performed using a supercomputing grant from the National Center of Atmospheric Research, which is sponsored by the National Science Foundation. Computational resources were supported in part by a grant of HPC time from the Department of Defense Major Shared Resource Centers.
BLASIUS Pressure Gradient Formulation
The sensitivity of mountain-wave solutions to the approximations made to remove acoustic waves from the equations of motion has been previously been demonstrated (e.g., Nance and Durran 1994; Nance 1997). The relatively rapid growth of wave amplitude with height in the BLASIUS simulations can be shown to be due to the formulation of the pressure gradient term in the anelastic equations used in this model. The BLASIUS equations can be written as
where u is the velocity, θ is the potential temperature, θ′ is the perturbation potential temperature, p′ is the perturbation pressure (to a reference, hydrostatically balanced profile), g is the gravitational acceleration, and and are the reference density and potential temperature profiles, respectively, both of which are functions of height only. In contrast, other commonly used forms of the anelastic equations involve writing the pressure gradient term as (e.g., Derbyshire et al. 1994). While the appearance of density within the derivative (rather than outside) is perhaps initially surprising, it arises from a more accurate treatment of the buoyancy term whereby the buoyancy force includes a contribution from the pressure perturbation. In the BLASIUS equation set this contribution is neglected. Note that in the frequently used Lipps and Helmer (1982) anelastic equations the gradient term in the momentum equation takes the form , where π′ is the perturbation to the Exner pressure and cp is the specific heat capacity at constant pressure. Again, unlike the BLASIUS equations, the reference-state variable appears within the derivative rather than outside it.
The consequences of adopting different forms of the pressure-gradient term in the anelastic momentum equation becomes apparent when the vertical structure equation for small-amplitude (linear) stationary mountain waves is derived. For simplicity we restrict attention to steady flow over a two-dimensional ridge and to a height-independent background flow . By splitting the velocity field into the sum of the background flow and small-amplitude perturbations, , substituting into Eqs. (A1)–(A3) and neglecting nonlinear terms, the following vertical structure equation can be obtained:
where k is the horizontal wavenumber and is the Fourier transform of the vertical velocity w′ defined as
Solutions to Eq. (A5) take the following form:
where A is a constant and . Thus, for a height-independent background flow, small-amplitude mountain-wave solutions to the BLASIUS equations will grow with height at a rate proportional to . Similarly it can be shown that when the pressure-gradient term in Eq. (A1) is written as , the vertical structure equation becomes
and the form of the second derivative in Eq. (A7) is different to that in Eq. (A4). Note that an identical equation to Eq. (A7) is also obtained for the commonly used Lipps and Helmer (1982) anelastic equations. Solutions to Eq. (A7) take the following form:
and is the density-scale height. The terms involving in Eq. (A9) are in general relatively small and result in only subtle changes to the vertical wavenumber m from that for the BLASIUS equations. However, the growth rate of the wave amplitude is significantly slower, depending as it does on rather than . We should therefore expect that in general mountain-wave solutions to the BLASIUS equations will grow more rapidly with height than those for anelastic equations where the pressure gradient takes the form or for the Lipps and Helmer equations.
Additional affiliation: Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, Colorado.
Current affiliation: NOAA/National Severe Storms Laboratory, Norman, Oklahoma.
This article is included in the Terrain-Induced Rotor Experiment (T-Rex) special collection.
COAMPS is a registered trademark of the Naval Research Laboratory.