The sensitivity of downslope wind forecasts to small changes in initial conditions is explored by using 70-member ensemble simulations of two prototypical windstorms observed during the Terrain-Induced Rotor Experiment (T-REX). The 10 weakest and 10 strongest ensemble members are composited and compared for each event.
In the first case, the 6-h ensemble-mean forecast shows a large-amplitude breaking mountain wave and severe downslope winds. Nevertheless, the forecasts are very sensitive to the initial conditions because the difference in the downslope wind speeds predicted by the strong- and weak-member composites grows to larger than 28 m s−1 over the 6-h forecast. The structure of the synoptic-scale flow one hour prior to the windstorm and during the windstorm is very similar in both the weak- and strong-member composites.
Wave breaking is not a significant factor in the second case, in which the strong winds are generated by a layer of high static stability flowing beneath a layer of weaker mid- and upper-tropospheric stability. In this case, the sensitivity to initial conditions is weaker but still significant. The difference in downslope wind speeds between the weak- and strong-member composites grows to 22 m s−1 over 12 h. During and one hour before the windstorm, the synoptic-scale flow exhibits appreciable differences between the strong- and weak-member composites. Although this case appears to be more predictable than the wave-breaking event, neither case suggests that much confidence should be placed in the intensity of downslope winds forecast 12 or more hours in advance.
Downslope windstorms are a frequently occurring phenomenon capable of adversely impacting commerce and communities in the lee of mountain barriers. Although high-resolution, cloud-resolving numerical weather prediction (NWP) models have shown an ability to realistically hindcast downslope windstorms (e.g., Colle and Mass 1998; Doyle and Shapiro 2000), the predictability of downslope winds remains largely unknown. Lorenz (1969) argued that, because of the rapid upscale propagation of uncertainties in the specification of initial conditions, the predictability of forecasts for mesoscale motions with spatial scales on the order of 10 km would be limited to time scales on the order of 1 h. This discouraging prospect has largely been supplanted by a more optimistic view based on experiences with high-resolution NWP models demonstrating that realistic mesoscale circulations can be generated during the forecast without having to specify mesoscale precursors to these circulations in the initial conditions. Anthes et al. (1985) argued that the predictability of many mesoscale phenomena can substantially exceed that suggested by Lorenz if they are sufficiently organized by the large-scale flow (e.g., fronts) or controlled by well-known external forcing such as orography, thermal contrasts, or other land-use characteristics.
Klemp and Lilly (1975) presented some of the earliest results suggesting that topography may enhance mesoscale predictability. They predicted downslope winds in Boulder, Colorado, using a linear two-dimensional mountain-wave model initialized with soundings taken from large-scale NWP model forecasts at locations upstream of the Colorado Rockies. Using soundings valid 3–5 h prior to the windstorm, they obtained a significant positive correlation between the predicted and observed wind speeds, at least for those events where gusts in excess of 25 m s−1 were measured. The results of Klemp and Lilly offer the possibility that downslope windstorms “are potentially predictable well in advance of their formation given accurate large-scale forecasts” (Anthes 1984); however, the question of what constitutes a sufficiently accurate large-scale forecast has remained largely unexplored.
A more recent attempt to forecast downslope winds produced more equivocal results. Nance and Coleman (2000) employed a two-dimensional nonlinear mesoscale model to forecast downslope winds using initial conditions generated by a synoptic-scale NWP model. Their approach, which was considerably more sophisticated than that used by Klemp and Lilly (1975) and which was based on synoptic-scale forecasts that should have been much better than those available in the 1970s, did have considerable success in providing alerts on days when windstorms were observed, but it also generated a large number of false warnings, suggestive of a difficulty in distinguishing between actual windstorms and null events.
A few previous studies have investigated aspects of downslope windstorm predictability using idealized two-dimensional numerical models. Doyle et al. (2000) compared simulations produced by 11 different numerical models initialized with identical upstream soundings representative of the conditions during the 11 January 1972 Boulder windstorm (Lilly and Zipser 1972; Lilly 1978). Although all of the models were able to produce downslope winds and a region of wave breaking in the lower stratosphere, the temporal evolution of the breaking was significantly different between the individual models. More recently, Doyle and Reynolds (2008) constructed an ensemble of two-dimensional simulations by adding random perturbations with magnitudes typical of radiosonde observational errors to the same Boulder sounding. For strongly nonlinear flow, the ensemble spread was characterized by a broad distribution with half of the members simulating a large-amplitude breaking wave in the stratosphere and the other half of the members simulating a trapped-wave response. Additionally, the downslope wind speed varied by over 25 m s−1 between the ensemble members.
The purpose of this study is to explore the predictability of downslope winds in a fully three-dimensional NWP model by examining the initial-condition sensitivity of two prototypical storms that occurred during the Terrain-Induced Rotor Experiment (T-REX; Grubišić et al. 2008). For each case, we conduct 70-member mesoscale ensemble simulations in which the initial-condition variability across the ensemble represents the expected level of uncertainty in the NWP analysis. We explore the time scale over which small initial differences between the ensemble members grow to large differences in the strength of the downslope winds. As we shall see, this time scale turns out to be far shorter than the one- or two-day period over which Anthes (1984) suggested accurate downslope wind forecasts could be made.
The two prototypical T-REX windstorms arise from different dynamical mechanisms. In the first case [intensive observation period (IOP) 6], a large-amplitude mountain wave breaks and strong downslope winds develop beneath the breaking region (Clark and Peltier 1977; Peltier and Clark 1979). Wave breaking can also lead to severe clear-air turbulence (CAT), which can be a significant threat to aviation (Lilly 1978; Nastrom and Fritts 1992). The mechanism that is active in the second case (IOP 13), which we will refer to as static-stability layering, occurs in cross-mountain flows with a strong stable layer near mountaintop level, below a much deeper layer of relatively weak upper-tropospheric stability (Durran 1986). For this class of windstorms, wave breaking is not an important factor, but severe downslope winds are still found along the lee slope.
The outline of this paper is as follows: the experimental setup, numerical model, and ensemble technique are described in section 2. The synoptic-scale evolution of the two events is given in section 3. The growth of ensemble variability for the two downslope wind events is described in section 4. Section 5 characterizes the challenges and limitations associated with downslope wind forecasting. Conclusions are presented in section 6.
2. Experimental design and setup
This section contains a description of the experimental design, numerical model, and ensemble method used in this study. A 70-member ensemble is used to characterize the initial-condition error growth for two strong downslope wind events which occurred during T-REX. T-REX took place from 1 March to 30 April 2006 over the southern Sierra Nevada and the adjacent Owens Valley (Fig. 1), and it sampled several strong mountain-wave and downslope wind events.
The first event we consider, IOP 6 (25–26 March 2006), generated a large-amplitude breaking mountain wave, along with vigorous turbulence and strong downslope winds. The second event, IOP 13 (16–17 April 2006), was characterized by upstream layering of the static-stability profile. Wave breaking was not apparent in the simulations for this case; however, strong downslope winds still developed.
a. Numerical model
Numerical simulations are computed with the Coupled Ocean–Atmosphere Mesoscale Prediction System (COAMPS; Hodur 1997), which solves a finite-difference approximation to the fully compressible, nonlinear, nonhydrostatic equations of motion on a terrain-following vertical coordinate system. Prognostic equations for the zonal, meridional, and vertical velocities (u, υ, and w), as well as potential temperature θ and Exner function π, are marched in time with a two-time-step integration (Klemp and Wilhelmson 1978). The finite-difference approximations to the space derivatives are second-order accurate, with the exception of the horizontal-advection terms, which are fourth-order accurate. Nonlinear instability is controlled in the model with a fourth-order diffusion operator applied to the prognostic variables at each time step. A full suite of physical processes are represented in the COAMPS simulations, including the parameterization of the boundary layer, radiative transfer, and moist processes (see Hodur 1997).
To resolve the multiscale aspect of downslope winds and mountain waves, three one-way interacting nests with 27-, 9-, and 3-km horizontal resolution are centered over the Sierra Nevada and Owens Valley (Fig. 1). The outermost 27-km mesh encompasses the western North American continent as well as a large portion of the eastern Pacific Ocean (Fig. 1a). An intermediate-sized 9-km mesh extends over the north–south extent of the Sierra Nevada range (Fig. 1b). The 3-km mesh is located over the highest portion of the Sierra Nevada and covers the region immediately upstream and downstream of the mountains, including the Owens Valley (Fig. 1c). Each nest is configured with 40 vertical levels, typical of real-time operational NWP models. The time step on the 3-km domain is 3.3 s and increases by a factor of 3 for each larger domain.
The lateral boundary conditions on the 27-km domain are specified from operational forecasts of the Naval Operational Global Atmospheric Prediction System (NOGAPS) and perturbed for each ensemble member using the fixed-covariance-perturbation method described in Torn et al. (2006). The perturbations are Gaussian with zero mean and a covariance that is balanced for synoptic-scale motions. The covariance relations are obtained from the Weather Research and Forecasting-Variational (WRF-VAR) data assimilation system (Barker et al. 2004). The 27-km domain was designed to be large enough so that these perturbations undergo several assimilation cycles before impacting the Sierra Nevada and Owens Valley. The lateral boundaries on the 9- and 3-km domains are updated every time step from the parent ensemble member. Topography on the lower boundary is specified from the global land one-kilometer base elevation (GLOBE) dataset and interpolated to the model grid points. A 25-point filter is applied to the interpolated topography field to remove the 2Δx signal. Upward propagating gravity waves are damped over the top 7 model levels by gradually increasing the numerical viscosity.
b. Ensemble method
An ensemble Kalman filter (EnKF) is used for data assimilation and to provide a set of 70 initial conditions for the numerical simulation of the downslope wind events. Under the assumption of Gaussian error statistics, the EnKF optimally combines a previous forecast with new observational data (Hamill 2006). The variability within the ensemble analysis represents the expected level of error, given uncertainty in the observations and in the previous forecast. The growth of ensemble variability over a given forecast period reflects the potential for error growth resulting from initial-condition uncertainty. A brief description of our implementation of the EnKF on the outer 27-km domain is given in the appendix. For detailed reviews of data assimilation using an EnKF, the reader is referred to Evensen (2003) and Hamill (2006).
For several reasons, it was beneficial to cycle the EnKF data assimilation system on the 27-km COAMPS domain for the duration of the T-REX field campaign (1 March–30 April). The two-month-long integration generated a reasonably large set of statistics from which the EnKF could be evaluated and properly calibrated. Additionally, because of the high computational cost of integrating the high-resolution ensemble, only active mountain-wave periods could be simulated with the 9- and 3-km-domain ensemble.
The ensemble initiation procedure for the 9- and 3-km domains for both the IOP 6 and IOP 13 simulations is shown schematically in Fig. 2. For the IOP 6 simulation, each ensemble member from the 27-km domain was linearly interpolated to the 9-km domain at 0600 UTC 24 March, and the 9-km ensemble was interpolated to the 3-km domain at 0000 UTC 25 March. From the interpolation time (when each nested domain is initialized) to 1800 UTC 26 March, each nest is independently cycled every 6 h with EnKF data assimilation. The predictability experiment begins at 1800 UTC 26 March by integrating the ensemble forward for 12 h without further data assimilation. In a similar manner, but for the IOP 13 simulation, the 9-km domain was initialized at 1200 UTC 15 April by interpolating data from the 27-km domain, and the 3-km domain was initialized at 0600 UTC 16 April by interpolating data from the 9-km domain. The set of initial conditions for 1800 UTC 16 April is integrated forward for 12 h without data assimilation. Interpolating the 27-km-domain ensemble to the finer meshes was a natural way to initialize the high-resolution ensemble forecasts. This procedure allows the resolved circulations and covariance relations sufficient time to spin up on each nest and ensures that the initial conditions on the nest contain horizontal scales appropriate for that grid spacing. Doyle and Smith (2003) have shown that mountain-wave simulations only require a few assimilation cycles to adequately represent the scales of motion resolvable by the numerical mesh. Similar nesting procedures for ensemble data assimilation have been performed in simulations of tropical cyclones (Zhang et al. 2009; Torn and Hakim 2009).
3. Synoptic-scale flow
The synoptic-scale flow pattern for the IOP 6 and IOP 13 events are described in this section. For IOP 6, the 500-hPa ensemble-mean analysis of geopotential height and wind speed is shown in Fig. 3. At 1800 UTC 25 March 2006 (Fig. 3a), a relatively sharp, negatively tilted trough is situated directly over the northern California coastline. A 45 m s−1 jet maximum is located at the base of the trough with strong southwesterly flow extending northeastwards to the northern Sierra Nevada. Six hours later, at 0000 UTC 26 March 2006 (Fig. 3b), the trough has progressed eastward such that the strongest winds are directly upstream of the Owens Valley and oriented nearly perpendicular to the barrier. The spatial extent of the high wind speeds has decreased and their maximum intensity has dropped below 45 m s−1. As will be seen, the ensemble mean downslope wind response is very strong at this time.
The synoptic-scale evolution for the IOP 13 case is shown in Fig. 4. At 0000 UTC 17 April 2006 (Fig. 4a), a broad, positively tilted trough is situated nearly 400-km west of the San Francisco Bay. A broad region of strong zonal flow, with wind speeds greater than 45 m s−1, extends from the base of the trough to the central and southern Sierra Nevada. Six hours later (0600 UTC 17 April 2006; Fig. 4b) the trough has intensified and progressed eastward to the upstream side of the Sierra Nevada. At this time, the flow is perpendicular to the Sierra Nevada crest; as will be seen in the next section, distinct static-stability layering is present throughout the troposphere, with strong stability in the lower troposphere and weaker stability farther aloft.
4. Downslope wind variability
In this section, the predictability of the wave-breaking event (IOP 6) and that of the layered event (IOP 13) is explored by considering the growth in the variability of the downslope wind response. To quantify the downslope wind intensity, an average zonal wind is computed on the 3-km domain within the white box depicted in Fig. 1c. To capture the low-level westerly momentum associated with the downslope winds, the metric box extends from ground level to 350-m AGL. This metric is referred to as the “Owens Valley metric” and is used to characterize the downslope wind response throughout the rest of this paper. The box covers the approximate horizontal extent of the T-REX observational network in which downslope winds were typically observed. For both cases, we consider forecasts valid at the time when the ensemble mean predicts a maximum downslope wind response; that is, 0000 UTC 26 March 2006 for the IOP 6 wave-breaking case and 0600 UTC 17 April 2006 for the IOP 13 layered case.
a. Ensemble distributions
To begin evaluating predictability characteristics of the two windstorm forecasts, we consider ensemble-derived distributions of the Owens Valley metric. Figures 5a,b show ensemble distributions of the Owens Valley metric for the 1800 UTC 25 March IOP 6 and 0000 UTC 17 April IOP 13 analyses, respectively. The distributions are computed by binning the Owens Valley metric from each ensemble member into 2.5 m s−1 wide bins and normalizing so that the area under the curve is unity. Equivalently, these can be thought of as ensemble-derived probability density functions (PDFs). These distributions show the variability of the lee-slope winds in the 70-member ensemble at the start of the forecast.
As a consequence of data assimilation, where observations systematically reduce the ensemble variance, the variability of the Owens Valley metric at the initial time is relatively small. For the 1800 UTC IOP 6 analysis (Fig. 5a) the majority of ensemble members lay between 15 and 22.5 m s−1, with a definitive peak in the 20 m s−1 bin. In contrast, the distribution in the IOP 13 analysis at 0000 UTC (Fig. 5b) is broader with the wind speeds at the initial time ranging between 5 and 17.5 m s−1 and a maximum in the 12.5 m s−1 bin. The initial distributions for both the 1800 UTC IOP 6 case and the 0000 UTC IOP 13 case are approximately normally distributed, which is consistent with the EnKF assumptions.
The PDF for the Owens Valley metric of the 6-h IOP 6 forecast (valid 0000 UTC 26 March) is shown in Fig. 5d. Compared to the initial time, the 6-h forecast contains a large level of uncertainty with a broad region of relatively uniform probabilities between 22.5 and 42.5 m s−1. Assuming that the observed response is sampled from the same distribution as the ensemble, there is approximately a 15% chance that the forecast winds will verify within any 5-m s−1 wide band over this range. In addition, a long tail extends toward weaker winds, suggesting a relatively high chance for false-positive forecasts. On the strong side, there is a sharp cutoff of the distribution for wind speeds greater than 42.5 m s−1; however, one very strong ensemble member is predicting downslope winds close to 50 m s−1.
The extent to which the spread in the forecast PDF might be due to small differences in the timing of the winds is explored in Fig. 6, which shows a PDF of the maximum values of the Owens Valley metric computed over a 6-h window centered on the 6-h IOP 6 forecast. The long tail extending from 7.5 to 20.0 m s−1 apparent in the PDF of instantaneous winds (Fig. 5d) is not present in this distribution. Nevertheless, there is still considerable spread with the predicted downslope wind response ranging from 20.0 to 50 m s−1. Even a downslope wind forecast verifying over a relatively long 6-h window would contain considerable uncertainty.
The 6-h IOP 13 forecast distribution (valid 0600 UTC 17 April) is shown in Fig. 5e. Although the wind intensities in the distribution range from nearly calm to a severe downslope windstorm with winds exceeding 25 m s−1, there is a distinct peak centered at the 17.5 m s−1 bin. The probabilities within this bin are nearly twice as large as the probabilities in the surrounding bins, suggesting more confidence can be placed in the 6-h forecast of the IOP 13 downslope wind event than in the 6-h forecast of the IOP 6 event.
The relatively narrow probability distribution for 6-h IOP 13 forecast suggests that it may be worthwhile to examine a 12-h forecast for the same event. The distributions for the analysis at 1800 UTC 16 April and the 12-h forecast valid at 0600 UTC 17 April are shown in Figs. 5c,f. As with the other two simulations, the analysis is characterized by a relatively narrow distribution associated with the data assimilation procedure; the wind speeds are contained within the 0 to 10 m s−1 range (Fig. 5c). However, over the course of the 12-h simulation, the wind speeds increase substantially in a number of the ensemble members and considerable uncertainty is associated with this downslope wind forecast (Fig. 5f). A relatively large section of the PDF is uniformly distributed, with approximately a 20% chance that the actual wind speeds will fall within any 5 m s−1 band between 2.5 and 17.5 m s−1. Compared to the 6-h forecast valid at the same time, a relatively long tail extending from 22.5 to 35 m s−1 indicates that a nontrivial chance of severe winds exists.
In the following, we will focus on the 6-h IOP 6 forecast initialized at 1800 UTC 25 March 2006 and the 12-h IOP 13 forecast initialized at 1800 UTC 16 April. These two cases will be referred to as the IOP 6 and the IOP 13 forecasts, respectively.
b. Strongest and weakest members
The predictability of the two events can be further examined by considering the range over which the forecast downslope winds vary within the ensemble. To this end, the members are ranked according to the forecast intensity of the Owens Valley metric, and the 10 strongest and 10 weakest ensemble members are grouped into two subsets. The shaded regions of the ensemble forecast distributions (Figs. 5d–f) show the fraction of the PDF containing the 10 strongest and 10 weakest members. Although these subsets are located on the tails of the distributions, they represent almost 30% of the total probability.
Averages over these subsets are computed to give composites of the strongest and weakest responses. These composites facilitate a simple comparison between the strong and weak cases on both the mesoscale and the synoptic scale. Furthermore, using 10 ensemble members to compute the composites ensures that our results are not unduly influenced by one or two unrepresentative outliers.
The Owens Valley metric at hour 6 (corresponding to 0000 UTC 26 March) was used to rank the 70 ensemble members for IOP 6. Figure 7 shows the evolution of the Owens Valley metric for each of the 10 strongest and 10 weakest ensemble members, as well as the mean for each 10-member subset. At hour 6, the mean for the strong subset is nearly 29 m s−1 stronger than the mean for the weak subset; however, the differences between the subset means decrease rapidly after hour 6. Looking at the individual members at hour 6, every strong member has wind speeds greater than 35 m s−1; whereas, the speed for almost every weak member is less than 15 m s−1. As will be discussed later, this large variability is associated with mountain-wave breaking above Sierra Nevada crest. At later times, the Owens valley metric varies dramatically among the individual members of both subsets, with wind speeds ranging from 5 to 35 m s−1.
Consistent with our focus on the 12-h forecast in the somewhat more predictable IOP 13 event, the weakest and strongest 10-member subsets are determined by the value of the Owens Valley metric at forecast hour 12 (0600 UTC 17 April). Figure 8 shows the individual ensemble members as well as the ensemble mean from the strong and weak subsets. Both the strong and weak subsets predict relatively intense downslope winds 3.5 h into the simulation that subsequently diminish by hour 4.5. Between 5 and 12 h, the differences between the two subsets steadily grows, with the strong-subset mean increasing to 26 m s−1 and the weak-subset mean decreasing to just 4 m s−1.
c. Downslope wind and mountain-wave response
A large-amplitude mountain wave generates the strong downslope winds in the Owens Valley. Here, the difference in the mountain-wave structure associated with the strongest and weakest ensemble members is compared along vertical cross sections depicted by the solid black lines in Fig. 1c. The orientation of the cross section is representative of the ensemble-mean midtropospheric synoptic-scale flow for the 6-h IOP 6 forecast (A–A’ cross section) and the 12-h IOP 13 forecast (B–B’ cross section). Both the A–A’ and B–B’ cross sections pass through the center of the Owens Valley metric box.
1) Wave breaking: IOP 6
Figures 9a,b show the zonal-wind component, as well as the turbulent kinetic energy (TKE) for the 6-h IOP 6 forecast along the A–A’ vertical cross section. The panels labeled “weak members” and “strong members” are obtained by averaging the fields over the weak and strong ensemble subsets. Consistent with the evolution of the Owens Valley metric shown in Fig. 7, the differences between the means of the strong and weak subsets have grown very large over the short 6-h forecast. For the strong subset (Fig. 9b), a tongue of high winds extends from the midtroposphere down the lee slope of the Sierra Nevada and into the Owens Valley. An extensive region of wave breaking is indicated by the strongly decelerated flow and large area of turbulent mixing between 8 and 12 km. The wave breaking is associated with the generation of the strong downslope flow on the lee slope (Clark and Peltier 1977; Peltier and Clark 1979). In contrast, the weak members (Fig. 9a) are characterized by high zonal-momentum air that does not extend below crest level. For this subset, the upper-level wave breaking is less extensive, weaker, and displaced vertically. A small region of decelerated flow, with zonal winds less than 10 m s−1 is apparent near 12 km MSL. Furthermore, the spatial extent of the upper-tropospheric turbulent mixing region is much smaller, with the TKE barely exceeding 10 m2 s−2.
The subset-mean vertical velocities w and potential temperatures θ for the 6-h forecast in IOP 6 are plotted along the same vertical cross section in Figs. 10a,b. The strong-subset mean (Fig. 10b) contains a large-amplitude mountain wave, as indicated by the large excursions of the isentropes and the couplet of intense vertical velocity on the lee side of the Sierra Nevada. For example, the 320-K isentrope is displaced downward nearly 4 km from its nominal upstream height of 8 km, and the maximum vertical velocity is nearly 14 m s−1 in the core of the updraft. Several of the individual ensemble members predict a much stronger mountain wave with vertical velocity magnitudes greater than 26 m s−1. Additionally, overturning is more evident for the individual members; however, the subset average smoothes this response. In contrast, the amplitude of the mountain wave for the weak-subset mean (Fig. 10a) is little more than half that of the strong subset. For the weak subset, the maximum vertical velocity is 8 m s−1 and the downward displacement of the 320-K isentrope is less than 2 km.
It is important to note that the difference between the strong and weak downslope wind response is due to major structural differences in the mountain wave, as opposed to small shifts in the leeward extent of some region of high surface winds. These large differences occur even though the upstream conditions are very similar. For example, the zonal momentum upstream of the Sierra Nevada crest is nearly indistinguishable between the strong and weak ensemble members (cf. the left-hand sides of Figs. 9a,b). Forward shear is apparent in both examples, with the zonal-wind increasing from 10 m s−1 near crest level to 40 m s−1 near the tropopause. The stability of the upstream profiles is also very similar between strong and weak members (Figs. 10a,b). A layer of strong crest-level stability is present and the undisturbed tropopause height is nearly identical for both subsets. Despite these similarities, the difference between the mountain waves and downslope wind forecast for the two subsets is considerable, suggesting very strong sensitivity to the model initial conditions and predictive time scales of less than the 6-h length of the forecast.
CAT commonly affects aircraft over regions of complex terrain (Nastrom and Fritts 1992). With finescale numerical models, it is possible to simulate CAT (Clark et al. 2000); however, the predictability of these events is not known. The large region of turbulent mixing and strong vertical velocities associated with the strong subset suggest that significant clear-air turbulence is forecast through a deep layer over the Owens Valley. In contrast, the weak vertical velocities and absence of TKE indicate that significant CAT is not forecast for the weak members. Thus, in this event, the predictability of CAT appears to be just as limited as the downslope wind predictability.
2) Stability layering: IOP 13
Consider now the strong- and weak-subset means for the IOP 13 case. Figures 9c,d shows the zonal wind, as well as the TKE for the 12-h IOP 13 forecast along the B–B’ vertical cross section. Evident in the strong-subset mean (Fig. 9d) is a region of high westerly momentum air extending down the lee slope of the Sierra Nevada and into the Owens Valley. The strength of the downslope flow exceeds 40 m s−1 high on the lee slope but decreases sharply toward the base of the Owens Valley. In contrast to the IOP 6 strong subset, mountain-wave breaking is not present in the upper troposphere along the B–B’ cross section. The upper-level winds above the Owens Valley are not reversed, but they are generally greater than 20 m s−1, and turbulent mixing, as indicated by TKE, is completely absent. Despite the lack of wave breaking, the difference between the strong and weak subsets is large. Strong zonal flow is only present well above crest level in the weak ensemble members (Fig. 9c). Additionally, the weak members exhibit easterly flow within the Owens Valley, which extends nearly halfway up the lee slope of the Sierra Nevada.
The θ and w fields in the mountain wave for the 12-h IOP 13 forecast are plotted along the B–B’ vertical cross section in Figs. 10c,d. In the strong-subset cases (Fig. 10d), the mountain wave generates trough-to-crest displacements of the 320-K isentrope of roughly 3 km and a vertical velocity maximum greater than 12 m s−1. Upstream of the Sierra Nevada, the static-stability profile is composed of two distinct layers within the troposphere: a layer of strong static stability below 6 km MSL and layer of weak static stability between 6 and 11 km MSL. The layered structure leads to an amplification of the shorter wavelengths through nonlinear processes (Durran 1986, 1992); in fact, the horizontal wavelength is roughly half that of the strong IOP 6 solution (Fig. 10b).
On the other hand, the weak-subset-mean mountain wave has much lower amplitude (Fig. 10c). The trough-to-crest displacements are reduced to roughly 1 km and the vertical velocity maximum is less than 8 m s−1. Upstream of the Sierra Nevada, a layered stability structure is again evident, but the vertical distribution of the Brunt–Väisälä frequency and wind speed is considerably different. Compared to the strong subset, the stability is weaker near crest level and stronger in the middle to upper troposphere and the wind speeds are greater in the upper troposphere. These differences in static stability and wind speed are likely responsible for the differences in the strength of the mountain wave and the downslope winds. The fact that such a large variation develops between the subsets over just a 12-h forecast suggests that windstorms produced by stability layering can also be quite sensitive to the model initial conditions.
5. Variability in the large-scale flow
In this section, the synoptic-scale variability between the strong- and weak-member subsets will be examined both at the time of the most intense downslope winds and one hour prior to that wind maximum.
a. Upstream soundings
A single sounding profile upstream of a mountain barrier is sometimes used for downslope wind prediction (e.g., Klemp and Lilly 1975; Nance and Coleman 2000). Here, we compare the upstream profiles of the cross-barrier component of the flow U, potential temperature θ, Brunt–Väisälä frequency N, and relative humidity RH between the strong and weak subsets of the IOP 6 and IOP 13 forecasts. For both cases, we consider model soundings computed on the 3-km domain one hour prior to the time in which the ensemble members were ranked. To determine appropriate locations for the upstream soundings, 1-h back trajectories were launched on the 3-km domains from the upstream edge of the Owens Valley metric box at a height of 5 km MSL, thereby ensuring that the profiles represent the atmospheric conditions over the crest at the time when the ensemble members are ranked.
1) IOP 6
Figure 11 shows soundings of U, θ, N, and RH from the 5-h forecast (valid 2300 UTC 25 March) for the IOP 6 simulation. With the exception of the low-level RH, the strong and weak subsets are strikingly similar. For example, the U profiles differ by less than 3 m s−1 through the depth of the troposphere, whereas the θ and N profiles are nearly indistinguishable. The soundings for each subset contains a stable layer centered at a height of 4.2 km, weaker static stability in the upper troposphere, and a tropopause height close to 9.5 km. A short-vertical-wavelength fluctuation in N is apparent in the 1.0–2.5-km layer of the weak-subset sounding. The influence of this small-scale feature on the intensity of the lee-side response is unknown, but it is expected that the deeper layer of strong static stability extending through the midtroposphere is much more important in regulating the development of strong downslope winds. It is hard to imagine that a forecaster could differentiate between the strong and weak responses based on these upstream profiles of U, θ, and N. Furthermore, assimilating radiosonde data at this location would do little to improve an NWP downslope wind forecast because, except for the humidity, the differences between the weak and strong members are within the error bounds associated with radiosonde observations.
In contrast to the wind speed and temperature, substantial humidity differences are apparent in the 1–3-km layer. Throughout this layer, the RH for the strong-subset mean is generally greater than 80% with a maximum value of approximately 95% at 1.5 km. In contrast, for the weak-subset mean, the RH is less than 80% in the 1–3-km layer with a minimum near 65% at 1.75 km. The importance of these differences in moisture will be discussed later.
2) IOP 13
For the IOP 13 case, the upstream profiles of U, θ, N, and RH are shown in Fig. 12. The distinct layering of static stability in the troposphere is evident; however, the structure is considerably different between the strong and weak subsets. On one hand, for the strong subset, a classic layered structure of the static stability is present with N exceeding 0.014 s−1 between the elevations of 2.5 and 5.2 km and lower stability farther aloft. This layering is very favorable for strong downslope winds on the lee side of the barrier (Durran 1986). On the other hand, the weak members are characterized by weaker stability in the 2.5–5.2-km layer and much stronger stability farther aloft. In addition to the differences in the static-stability profiles, the cross-barrier flow above 6-km is nearly 10 m s−1 faster in the weak-member subset. In contrast to the other variables, the RH profile is very similar between the strong- and weak-member subsets, with significant drying above 3 km in both profiles. Thus, unlike IOP 6, signatures potentially capable of distinguishing between the strong and weak subsets are available at this short 1-h lead time, suggesting a somewhat longer predictive time scale for IOP 13.
3) Moist processes
The large difference in low-level moisture presented in Fig. 11d suggests the possibility that diabatic effects associated with moist processes are modulating the strength of the windstorm. Jiang and Doyle (2009) argued that low-level moisture during IOP 6 strengthened the mountain-wave response by reducing the below-mountaintop stability and decreasing the depth of the blocked layer. To asses the role of moisture and latent heating on the error growth of the Owens Valley metric a “dry” version of the IOP 6 ensemble simulation is conducted in which diabatic and radiative effects associated cloud microphysical processes are eliminated. In all other respects, including the set of initial and boundary-conditions perturbation making up the ensemble, the new set of simulations is identical to that of the original “control” IOP 6 experiment.
We first evaluate the downslope wind response for the same set of ensemble members, which comprised the strong- and weak-member subsets in the control simulation. Consistent with Jiang and Doyle (2009), the strong-member-subset mean weakens by 12.5 m s−1 at hour 6. On the other hand, the dry weak-member-subset mean strengthens by 3.6 m s−1. As a consequence, the difference between the means of the strong and weak subsets is reduced from 29 m s−1 in the moist case to 13 m s−1 in the dry case. Although the difference in wind speed remains substantial, it does appear that the low-level humidity plays a significant role in differentiating between the strong and weak cases in IOP 6.
The generality of this result is nevertheless unknown. The strong- and weak-member subsets in the IOP 13 case were not associated with large differences in upstream humidity. Furthermore, the difference between the strong- and weak-member-subset means in a dry simulation for the IOP 13 case only changed by 4.7 m s−1. In addition, if the ensemble members in the dry IOP 6 simulation are reranked to correspond to the 10 weakest and 10 strongest members in the dry case (instead of retaining the rankings from the moist simulation), the strong- and weak-member-subset means differ by 25 m s−1, which is only 4 m s−1 less than the difference obtained in the original control simulation. However, as shown in Fig. 13, the 1-h upstream soundings for the strong- and weak-member subsets are once again extremely similar, suggesting that, even in the absence of moisture, downslope winds generated by wave breaking can be extremely sensitive to small variations in the upstreamflow.
b. Synoptic-scale structure
Although the use of a single upstream sounding to forecast downslope winds may be attractive due to its simplicity, a more complete representation of the synoptic-scale data is likely to be more informative. Here we examine the differences in the synoptic-scale conditions associated with the weak and strong ensemble members by comparing subset-means on the 9-km-resolution domain.
As noted in section 3, both IOP 6 and IOP 13 occurred as progressive upper-level troughs crossed the Sierra Nevada. Figures 14a,b show the subset means of the 500-hPa wind speed V and geopotential heights from the 6-h forecast (valid 0000 UTC 26 March), the time of strongest downslope winds. The similarities between the strong and weak members is surprising. In both subsets, the trough is located directly over the Sierra Nevada crest and a wind speed maximum is just north of the Owens Valley. The difference in V is also evident in the vertical cross sections shown in Figs. 15a,b, which are taken along the line C–C’ indicated in Fig. 1b. Note that although the winds in the strong subset are slightly more intense north of the Owens Valley metric box (indicated by the dashed vertical lines), directly upstream of the metric box the difference in velocity is less than 2.5 m s−1. In contrast to the wind speed, there is almost no discernible difference in the θ field between the weak and strong subsets. Averaged over the entire cross section shown in Fig. 15, the RMS difference in V between the strong- and weak-subset means is 2.1 m s−1; for θ the RMS difference is 0.8 K.
Given the large variability in the simulated downslope wind response, it is remarkable that the differences in the synoptic-scale flow are so small. These differences are within the error bounds typically associated with radiosonde observations, further suggesting that differentiating between the strong and weak storms for this wave-breaking event would be nearly impossible for a human forecaster or, as seen in these simulations, a deterministic numerical model.
Turning now to the IOP 13 case, Figs. 14c,d show the 500-hPa wind speed and geopotential heights from the 12-h forecast (valid 0600 UTC 17 April) for the strong- and weak-subset means. The differences between the strong and weak members are clearly larger than in IOP 6. For example, compared to the weak-subset mean, the jet core in the strong-subset mean is displaced approximately 150 km northward and has a slightly greater horizontal extent. Relatively large differences between the strong and weak subsets are also apparent in the vertical cross sections of V and θ shown in Figs. 15c,d (which are again taken along the line C–C’). The RMS difference in V in the plane of the cross section is 7.5 m s−1, with the majority of the increase occurring north of the Owens Valley metric box. However, consistent with the upstream soundings shown in Fig. 12a, the winds directly upstream of the Owens Valley metric box are actually weaker in the strong subset. The RMS difference in θ is 4.1 K, which, as will be discussed later, is associated with vertical displacements in the frontal inversion.
The difference in the position of the 500-hPa trough between the strong and weak subsets is not a “timing error” because of differences in the propagation speed of the trough axis. The eastward propagation speed and location of the trough axis are nearly identical in both subsets. Instead, the differences arise from the meridional displacement of the jet stream and its associated upper-level front. In the strong subset, the northern location of the jet situates the sloping surface of the upper-level front near crest level directly upstream of the Owens Valley metric box. This configuration of high static stability in the lower troposphere and weaker static stability aloft provides a favorable configuration for the development of severe downslope winds. In contrast, the more southern location of the jet in the weak subset results in a southward displacement of the upper-level front. As a consequence, the high static-stability air associated with the frontal layer is farther aloft upstream of the Owens Valley metric box, resulting in an unfavorable layering structure of strong downslope winds.
c. Lateral boundary condition uncertainty
To evaluate the roll of lateral boundary condition uncertainty on the growth of the Owens Valley metric variability, two “deterministic boundary” experiments are performed in which the NOGAPS deterministic forecast is used to specify identical lateral boundary conditions for each ensemble member. We compare the mean of the Owens Valley metric for the same set of members that made up the strong- and weak-member subsets in the control experiment. In general, the boundary perturbations had only a small impact on the error growth of the Owens Valley metric. For example, in the IOP 6 deterministic boundary simulation, the difference between the strong- and weak-subset mean is 25.8 m s−1, compared to a difference of 28.3 m s−1 for the control ensemble. The IOP 13 simulations are just as insensitive to the boundaries perturbations. For this case, the difference between the strong- and weak-subset means for the same set of ensemble members as the control simulation is 24.3 m s−1, compared to a difference of 24.5 m s−1 for the control ensemble. Clearly, the lateral boundary perturbations do not significantly contribute to the ensemble variability with the limited predictability of these downslope windstorms being fully attributable to initial-condition uncertainties. Furthermore, inspection of the low-level zonal-wind ensemble variance on the 3-km domain (not shown) indicates that the large errors are primarily anchored to the lee side of the Sierra Nevada. This suggests that the rapid error growth of the Owens Valley metric occurs entirely over the Owens Valley and that large errors are not simply being advected through the 9-km domain lateral boundaries.
6. Summary and conclusions
We have attempted to systematically document the sensitivity of downslope wind forecasts to initial conditions in a three-dimensional NWP mesoscale model. An ensemble of 70 different initial conditions is generated for each of two prototypical downslope wind events from the T-REX special observing periods: IOP 6 and IOP 13. Consistent with the available data, most of the simulations of IOP 6 show a large-amplitude mountain wave with upper-level tropospheric wave breaking and severe downslope winds. In contrast, wave breaking was not present for the most of the IOP 13 simulations; consistent with the observations, the downslope winds were significant but weaker than in IOP 6. The strong winds in IOP 13 were generated by a layer of high static stability flowing beneath a mid- and upper-tropospheric layer of low stability.
In both cases, the individual ensemble members were ranked according to the forecast intensity of the near-surface winds in a region along the lee slope of the Sierras (the Owens Valley metric box), and the 10 strongest and 10 weakest ensemble members were grouped into separate subsets. For the wave-breaking simulations (IOP 6), initial-condition errors grow rapidly, leading to large variability in the downslope wind forecast. The differences between the means of the strong- and weak-member subsets grew from less than 5 m s−1 at forecast hour 3 to 28 m s−1 at forecast hour 6. These differences at forecast hour 6 represent the contrast between a very severe 41 m s−1 downslope wind and a mild 13 m s−1 event.
For the case with layered static stability but no wave breaking (IOP 13), initial-condition errors grow at a somewhat slower pace. For the 6-h forecast (valid at 0600 UTC 17 April), the difference between the strong- and weak-subset mean grows to 15 m s−1. Additionally, the ensemble-derived PDF of wind speed shows at least some degree of clustering about a mode of 17.5 m s−1. On the other hand, in the 12-h forecast (valid at the same time), the difference in the subset mean grows to 22 m s−1 over a period of roughly 8 h, and the PDF of the wind speed is much broader. The 26 m s−1 winds in the strong-subset mean constitute a moderate windstorm, whereas the 4 m s−1 winds in the weak subset are a nonevent.
Upstream soundings from the mean of the strong- and weak-member subsets were examined for both cases just one hour before the time of the maximum winds. In IOP 6, the wave-breaking case, the differences in the cross-barrier wind speed, potential temperature, and Brunt–Väisälä frequency for the two subsets were generally less than radiosonde observational errors. Dynamically significant differences were present in the low-level humidity field; however, even after accounting for the influence of moisture, the strength of the downslope winds remained very sensitive to small variations in the upstream profile. Such sensitivity suggests that deterministic operational forecasts may have difficulty accurately predicting the strength of downslope winds and CAT associated with mountain-wave breaking. This limited 1-h predictive time scale is consistent with the two-dimensional results for the wave-breaking regime obtained by Doyle and Reynolds (2008).
For the case with strong low-level static stability (IOP 13), in which wave breaking did not play a major role, the predictability time scale appears to be somewhat longer. Upstream soundings from 1 h prior to the time of maximum wind show clear differences between the mean profiles for the strong- and weak-member subsets. For example, a 2-km-deep layer of strong static stability is present directly above crest level in the strong subset, whereas the crest-level static stability is considerably lower in the weak subset. Downslope winds that develop in cases like IOP 13 may have a significant degree of predictability at 6-h lead times, but they may have much less at lead times of 12 h.
Some of the difficulty in distinguishing actual events from null cases encountered by Nance and Coleman (2000) when trying to forecast downslope winds with a mesoscale model may be due to the large initial-condition sensitivities documented in this paper. Although the rapid growth of uncertainty in our ensemble forecasts of the IOP 6 and IOP 13 events may not occur in some downslope windstorms, these results nevertheless call in to question the idea that the predictability time scale for downslope windstorms might be similar to that of the synoptic-scale flow. Particularly in the case of windstorms generated by wave breaking, the actual strength of the downslope wind may depend on small mesoscale circulations that represent only a modest component of the total synoptic-scale signal. Such mesoscale perturbations may be predictable over time scales that are no longer than those originally suggested by Lorenz (1969).
Finally, we note that we have focused on one fundamental aspect of predictability: the growth rate of small initial uncertainties. Other important considerations, such as the time required before a downslope wind forecast shows no more skill than persistence or climatology, are left for future research. For example, although in many respects the downslope winds in IOP 6 are very hard to predict, almost all 70 ensemble members produce some type of windstorm in the 6-h forecast. Given that high winds did actually occur, the 6-h ensemble forecast certainly exhibits significant skill relative to climatology. The skillfulness of the forecast may, however, be less clear cut if it is compared with the more limited climatology of lee-slope winds that occur when the cross-mountain flow at 700 hPa exceeds some modest threshold, such as 10 m s−1.
We thank Ryan Torn and Greg Hakim for the use of their EnKF code. Comments from two anonymous reviewers helped improve the manuscript. This research was supported by Office of Naval Research Grant N-00014-06-1-0827. Computational resources were supported in part by a grant of HPC time from the Department of Defense Major Shared Resource Center at Wright Patterson Air Force Base, OH. COAMPS is a registered trademark of the Naval Research Laboratory.
Implementation of the EnKF
Following Whitaker and Hamill (2002), a square root formulation of the EnKF is used. Because of computational constraints, the size of the ensemble was limited to 70 members; however, several studies have indicated that O(100) members are sufficient to perform data assimilation on the synoptic scale (Whitaker et al. 2004; Mitchell and Houtekamer 2002; Dirren et al. 2007; Torn and Hakim 2008). The 70-member ensemble was generated on the 27-km outer domain for the entire period of 1 March–30 April 2006. Observations from radiosonde data, Automated Surface Observing Systems (ASOS) data, Aircraft Communications Addressing and Reporting System (ACARS) data, as well as cloud-drift winds were assimilated into an ensemble of first guesses every 6 h. No attempt was made to assimilate special observations from the T-REX dataset on the 27- or 9-km domains; however, 10-m wind data were assimilated for three mesonet stations spread across the Owens Valley from the Desert Research Institute mesonetwork. The potential for large model errors associated with poorly resolved gravity waves (e.g., Reinecke and Durran 2009) limited the use of the 10-m wind data to the 3-km domain. The set of potential observations was limited by only considering data within a 2-h window centered on the assimilation time. Furthermore, cloud track winds were thinned by averaging all observations within a 1° and 25-hPa radius; ACARS data were thinned by averaging all data within a one-grid-point and 25-hPa radius; and surface observations were not used if the discrepancy between the station elevation and the model elevation was greater than 500 m. The initial 27-km-domain ensemble was constructed by interpolating a 36-h NOGAPS forecast, valid at 1200 UTC 1 March 2006, to the 27-km COAMPS domain and adding balanced random perturbations with zero mean and derived from the WRF-VAR data assimilation system (Barker et al. 2004). Following Torn and Hakim (2008), the initial perturbations were scaled by a factor of 1.75 so that the initial ensemble variance is slightly larger than the RMS error between the ensemble mean and radiosonde observations.
As a consequence of using a relatively small ensemble to sample the background-error statistics, the covariance estimates can generate spurious long-distance relationships (Houtekamer and Mitchell 2001; Hamill et al. 2001) as well as underestimate the magnitude of the covariance relations (Whitaker and Hamill 2002). To boost the rank of the ensemble, two standard procedures are used: covariance inflation and covariance localization. Long-distance correlations are restricted by applying a broad isotropic localization function [Gaspari and Cohn 1999, their Eq. (4.10)], which is unity at the observation location and decreases monotonically to zero at 3000 km from the observation. The localization radius was decreased by a factor of 3 for each nested domain. Covariance relaxation is used on the 27-km domain to boost the covariance magnitude (Zhang et al. 2004). The background- and analysis-error covariances are linearly averaged with weights of 0.775 and 0.225, respectively. These values were determined experimentally to ensure that the ensemble variance, averaged over the two-month period, was approximately equal to the RMS difference between the ensemble mean and radiosonde observations. Furthermore, to ensure that the ensemble was neither under- or overdispersive, rank histograms (Wilks 2006; Hamill 2001) were used to monitor the two-month-long integration (not shown). Because of the rapid growth of variability, as well as the short period of use, covariance inflation was not used for the 9- and 3-km ensembles.
* Current affiliation: Naval Research Laboratory, Monterey, California.
Corresponding author address: Patrick A. Reinecke, Naval Research Laboratory, 7 Grace Hopper Avenue, Monterey, CA 93943-5502. Email: email@example.com
This article included in the Terrain-Induced Rotor Experiment (T-Rex) special collection.