17.6

AN INVESTIGATION OF THE DIURNAL VARIABILITY OF THE

CENTRAL COLORADO DOWNSLOPE WINDSTORM

Daniel B. Weber *

Center for Analysis and Prediction of Storms

University of Oklahoma, Norman, Oklahoma

Strong westerly winds in the lee of the Central Colorado Rockies termed "downslope winds" have been observed for many years. This phenomenon also occurs near other mountain ranges such as the Andes in South America and the Alps in Europe and can be attributed to gravity wave-related processes (Long 1953, Smith, 1980, Lilly, 1983, Durran 1986a). Observational studies (Julian and Julian, 1969, Brinkmann, 1974, and Whiteman and Whiteman, 1974), using data collected in the Boulder Colorado area, suggest seasonal and diurnal variations in windstorm strength and occurrence. As depicted in Figure 1, a minima exists during the summer months of June, July, and August and a frequency maximum occurs during January. The annual peak in windstorm observations during January is associated with a minimum in solar radiation and strong cross-mountain flow. These studies indicate, as Figure 2 illustrates, that 70% of the reported windstorms occurred during nighttime hours. The data show a four-fold difference between the frequency maximum at night and the daytime minimum.

The meteorological literature contains hundreds of papers on gravity wave theory and hydraulic flow and their application to flow over a mountain. However, there are relatively few papers in the literature addressing the effects of a well mixed boundary layer on mountain wave flow (Raymond, 1972, Reisner and Smolarkiewicz, 1994). The purpose of this work is to investigate the observed windstorm tendencies via an idealized analytical solution and numerical simulations.

The simplified two-layer linear approach is chosen for two reasons. Linear theory captures the basic gravity wave structure and the two-layer configuration allows for the introduction of a neutral boundary layer. The resulting solution becomes a function of the mixed layer depth. ____________________________________________

Figure 1. Monthly windstorm frequency distribution for Boulder, CO. The data source is Julian and Julian (1969).

Figure 2. Hourly frequency distribution for windstorms observed in Boulder, CO during the period 1869-1972. The data is taken from Whiteman and Whiteman (1974).

In previous work involving surface heating, the heat source is located directly over the mountain and is not a part of the upstream condition. Therefore, the amount of time the parcel spends over the mountain is small compared to the total trajectory time. In the real atmosphere surface heating occurs far upstream of the mountain as well as near the mountain. On length scales of the mountain width, the top of the boundary layer is nearly horizontal.

A different approach is taken here in regards to the linear analysis. Steady state linear theory is applied from the standpoint that the convective motions, associated with the process of heat redistribution in the mixed layer are small compared to strong mountain wave circulations, are neglected. From this perspective, the solutions in the overriding stable layer can be easily solved in terms of the mixed layer depth and the horizontal and vertical scales.

During the first part of the diurnal heating cycle, the steady state assumption is not defensible as the mixed layer height is strongly time-dependant. The steady state assumption is valid in the late afternoon when the mixed layer height grows slowly and the heat is distributed over a large vertical extent (~1km). A time dependent solution in terms of the mixed layer height is not investigated here but is possible through the application of similarity theory (Garratt, 1994).

A linearized two-dimensional Boussinesq equation set in terms of , and is used in this analysis. Little generality is lost from the application of the Boussinesq set of equations, as the effects of decreasing density with height are well known (Gutman, 1991). For the analysis given here, the base state wind is constant () with height and equal in both layers. The stratification is defined as and is set to zero in the lower layer () and to a positive constant () in the upper layer. The steady state version of the Boussinesq equation set is:

, (1)

, (2)

, (3)

. (4)

Figure 3 sets up the problem graphically.

Layer 2

Layer 1

Figure 3. Graphical depiction of the two-layer linear problem. In the analysis, the mountain is chosen to be a cosine function.

Equations (1)-(4) combine to give a single equation in perturbation vertical velocity :

, (5)

with solutions of the form:

, (6)

. (7)

The wave numbers in the two layers are:

, .

The wave number in layer 2 is a function of the stability, base state wind, and horizontal wave number. Equation (6) is valid for real. If becomes imaginary then the solutions no longer admits gravity waves modes and follows (7). The wave number in the bottom layer is equivalent to the horizontal wave number. The complex coefficients A, B, C, and D are determined from application of the boundary and matching conditions. The upper boundary condition requires energy to propagate out of the domain. Following Eliassen and Palm (1960) this condition is met when . The overbar represents the horizontal average. The matching conditions require the displacement and pressure at the layer interface height (z=0) to be equivalent.

In addition, a discontinuity in density is introduced into the solution and represents an inversion placed at the interface between the two layers. An inversion is commonly observed at the top of the boundary layer. The application of the interface pressure condition follows that of Klemp and Lilly (1975) and can be obtained by integrating the hydrostatic relation. These conditions are presented in terms of the vertical velocity:

, at z = 0

, at z = 0

and z=0 is defined at the layer interface. The second matching condition requires continuity of the vertical derivative of the vertical velocity and is equivalent to matching the horizontal pressure gradient term in the two layers. The term represents the change in potential temperature across the inversion. Physically, the inversion represents external gravity waves along the layer interface. For the non-inversion case this term vanishes. The bottom boundary condition is linearized using:

= .

The terrain is defined by:

,

where, is the horizontal wave number and is the mountain height. Solving for the complex coefficients and incorporating them into (6) and (7) and taking the real parts gives:

, (9)

where the constants are given by:

, , ,

, , where .

Results from the two-layer solution are presented in terms of normalized surface wave drag. For a linear hydrostatic Boussinesq system, the steady state surface wave drag is equivalent to the vertical flux of horizontal momentum (Eliassen and Palm, 1960):

,

The vertical flux of horizontal momentum can be shown to be:

.

The sensitivity of the surface wave drag to the mixed layer depth, horizontal wavelength, and for the hydrostatic case the inversion strength, is illustrated by the curves in Figure 4. For this analysis, the base state wind is U= 20m/s and the static stability in the upper layer is N= 0.01. The flux is normalized by the H= 0 case. The steady state momentum flux in the upper layer is significantly affected by changes in the thickness of the mixed layer and to a lesser extent by the horizontal wave number. The results reveal a nearly 80% reduction of the wave activity in layer 2 for a mixed layer depth of 3km. The wave number dependence is small, with only a 6% decrease in the remaining wave activity for a 80km wavelength mountain case versus the 160km wide mountain flow for a 3km deep neutral layer. Interestingly, the plot also shows that an inversion acts to offset the reduction of mountain wave activity due to the development of a neutral layer. For a 1.5K inversion, a mixed layer of 0.5km depth is required to remove the inversion layer mountain wave enhancement. A 10K inversion, although not likely to be observed during windstorms, greatly enhances the mountain wave activity. The enhancement is likely due to surface wave effects, but a detailed study has not been performed to confirm this. This analysis can be extended to systems with more than two layers.

Figure 4. Plot of the linear analytical steady state surface wave drag curves as a function of mixed layer depth, inversion strength, and horizontal wavelength (80,160 km). The values are normalized by the H=0 steady state values. Plotted points represent the normalized surface drag at the conclusion of the heating cycle for the finite amplitude and 7 km heated critical layer simulations (from Sections 4.1 and 4.2).

A new three-dimensional model (ARPI3D, with a two-dimensional option) was used to simulate dry heated mountain wave flows. The model combines attributes from Durran and Klemp (1983) and ARPS (1995). It includes a simple coordinate transformation, a vertically-implicit solving technique, and a linearized upper radiation condition between pressure and vertical velocity. In addition, the model solves the more efficient non-dimensional Exner function instead of pressure. The primed quantities represent the perturbations from the steady base state (overbar). The equation set is:

with,

Normally, mesoscale models make a number of approximations that focus primarily on the linearization of the meteorologically insignificant acoustic modes. In this application, the non-linear pressure gradient terms were found to contribute significantly in strongly forced mountain wave flows. Perturbations on the order of 30% of the base state potential temperature were observed in the breaking wave regions in the stratosphere. The non-linear pressure gradient term reduces the time step slightly but improves the numerical solution. The non-linear divergence and adiabatic terms are included to close the model energetics. Little improvement is noted when the coefficient of the divergence term is at full strength.

The above equation set is transformed to a terrain following coordinate system according to Durran and Klemp (1983). The base state variables are horizontally homogeneous and hydrostatically balanced. Sub-grid scale motions are parameterized following the base model of Sullivan et. al. (1994). This approach utilizes a 1.5 order closure scheme that predicts turbulent kinetic energy for use in the determination of the eddy viscosity. During unstable conditions, the mixing length in the planetary boundary layer (PBL) is computed using the method described by Sun and Chang (1986). The surface heat flux is estimated with a sine function and the surface friction computed following Miller and Durran (1991).

The lateral boundary conditions for the normal velocity components are variants of that proposed by Klemp and Lilly (1979) in which the phase speed is estimated by the linear hydrostatic value. The top boundary condition follows the analytical approach of Klemp and Durran (1983) and allows upward propagating gravity wave energy to pass vertically out of the top boundary. The model equations are computed using spatially centered energy conserving finite differences on the Arakawa (1966) staggered C-grid. The time integration follows the split time step approach of Klemp and Wilhelmson (1978). The horizontal and vertical advection terms are differenced using an energy conserving scheme first proposed by Arakawa (1966) and later modified to the fourth order equivalent advective form by Xue and Lin (1991).

The numerical model was applied in two- dimensions to a series of non-linearly forced mountain wave flows and to an observed downslope windstorm. Each simulation was advanced to steady state conditions and the surface heating cycle initiated. Comparisons are made between the non-heated control run, their heated counterparts, and with linear theory. The mountain shape is estimated by a "Witch of Agnesi" profile.

Non-hydrostatic effects can be measured by the ratio of the horizontal to vertical wavelengths, where values approaching unity indicate scale equivalence and signify substantial non-hydrostatic forcing. For this test the ratio is:

,

where N=0.01956 is the static stability, =10 km the quarter wavelength parameter, U=20 m/s is the base state wind, and the mountain height h=300m . The non-linear effects are measured by the gravity wave strength (Nh) normalized by the base state wind.

.

The Figure 5 shows that the mountain wave strength decreases, measured by the vertical flux of horizontal momentum, as the depth of the mixed layer increases. The normalized surface wave drag computed at the end of the heating cycle is plotted in Figure 4 at the H=0.75km depth. The estimate from linear theory matches the numerical results.

Figure 5. Vertical profile of the vertical flux of horizontal momentum at Ut/a = 30, 60, 90, 120, and 150 for the heated mountain (a=10km) test. Profiles are normalized by the linear hydrostatic Boussinesq value. One vertical wavelength is 6.4 km.

Downslope windstorms frequently display elevated regions of flow reversals and enhanced turbulence with weak or neutral stability. These regions are commonly referred to as critical layers and exist when the phase velocity of the wave equals that of the transport medium. For the case of airflow over mountains this occurs when the cross-mountain wind speed is reduced to zero. A number of simulations were performed and the results compared to the control run, Durran (1986b), and to linear theory. In these simulations, a mean state critical layer was introduced at a height of 7km. The base state wind is constant (20m/s) from the ground to 5km, and then linearly reduced to zero at a height of 7km. Above 7km the wind speed is zero. The static stability is N=0.01047 and the inverse Froude number Nh/U = 0.314, 0.39, and 0.471 for the h=600m, 750m, and 900m mountain height tests, respectively.

As Figure 6 displays, the high drag states for the h= 750m and 900m tests remain at the conclusion of the heating cycle. These results compare favorably with Durran (1986b, Table 1). The application of parameterized heating reduces the steady state response but did not eliminate the high drag state. The reduction in wave drag from the control runs due to heating are plotted in Figure 4. The points representing the h=750m and 900m tests lie above the linear theory reduction curves. Linear theory over-predicts the decrease in wave activity by as much as 85%. The 7km h=750m 2 km deep neutral layer test requires nearly 20 hours to achieve an elevated drag state. The neutral layer test verifies the heated h= 750 m simulation but also exposes a slow yet significant growth mode.

Figure 6. Time series of computed surface wave drag for all tests with a mean state critical layer at 7 km. The solid lines represent the control runs and the dashed lines the heated simulations. The heating curve is provided at the bottom of the plot with a maximum heating rate of 200 and a minimum rate of -40(h=750m case only). The vertical lines depict the approximate mixed layer depth.

These tests include observations taken upstream of the Front Range and a realistic two-dimensional mountain profile. The model is initialized with the 2305 UTC atmospheric sounding collected from Craig, Colorado. As before, two runs were conducted, a control and a heated test. The maximum heat flux value was estimated from data obtained from the BOREAS experiment (Black, 1996).

The mountain induced surface wave drag and maximum surface wind speed, as a function of time, are presented in Figures 7 and 8, respectively. Both figures indicate a reduction in the windstorm strength. The wave drag computed from the 2305 UTC heated run (dashed line) is nearly 20% lower than the control run (solid line) at the conclusion of the heating cycle. The maximum surface wind speed is reduced by 15-20% from the control run and the depth of the mixed layer reached 1.5 km. The predicted maximum surface wind speed and observed wind measurements are not directly comparable, since the model did not include the change in the upstream conditions with time. But the simulated maximum surface wind speed is similar in magnitude to those measured near the end of the observed storm (not shown). Figure 9 compares the horizontal wind speed at the conclusion of the simulations.

Figure 7. Graph of the surface wave drag for the January 9, 1989 Boulder windstorm simulation for the control (solid line) and heated (dashed line) runs. The heating curve is provided at the bottom of the plot with a maximum heating rate of 200. Vertical lines indicate the approximate depth of the mixed layer.

Figure 8. Plot of the maximum surface wind speeds for the Boulder January 9, 1989 windstorm simulation for the control (solid line) and heated (dashed line) runs. The heating curve is provided at the bottom of the plot with a maximum heating rate of 200. Vertical lines indicate the approximate depth of the mixed layer.

(a) (b)

Figure 9. Horizontal velocity plots for the January 9, 1989 Boulder windstorm simulation at t = 70,000 seconds for (a) control and (b) heated tests. Area shown is the lower 6 km in the vicinity of the mountain peak and city of Boulder. The contour interval is 2.0 m/s.

In every test, parameterized surface heating reduces the established mountain wave activity, as measured by the surface wave drag, maximum horizontal velocity at the surface, and the vertical flux of horizontal momentum. Reductions on the order of 30% for the wave drag and 10-15% for the maximum surface winds were common in the heated simulations. The response to surface heating is found to be a function of the mixed layer depth, with deeper surface layers forcing larger reductions in the wave activity aloft. The largest differences between linear theory and the numerical results were noted in the critical layer tests where linear theory overestimated the actual wave reduction by 10% to 85%. In tests displaying wave-breaking characteristics the surface drag is diminished but a high drag state remains. The transition from the steady high drag state to a low drag state did not occur when influenced by a moderate amount of surface heating. The tests performed here applied surface heating at a rate similar in magnitude to that observed in Central Canada in February, consistent with the wintertime peak in windstorm observations. In all of the high drag simulations, the flow in the boundary layer is dominated by the upper layer wave response. This is expected since the depth of the mixed layer was rarely greater than 1/4 of a vertical wavelength. The surface wind speed did not undergo such a large decrease and is likely due to two factors. The drag is a quadratic quantity in perturbation variables and will respond to changes in the flow more rapidly. Secondly, potential flow theory requires an increase in the flow on approach to an obstacle.

These results support the Central Colorado observational windstorm record and point to a general cause in the windstorm strength and occurrence pattern.

I would like to thank my advisor Doug Lilly for his guidance and suggestions during this research and to CAPS, the University of Oklahoma, and the National Science Foundation for supporting this research.

Arakawa, A., 1966: Computational design for long term integration of the equations of motion: Two-dimensional incompressible flow. Part I. J. Comput. Phys. 1, 119-143.

ARPS, 1995: ARPS Users Guide Version 4.0, Center for Analysis and Prediction of Storms, University of Oklahoma, 380pp.

Brinkmann, W. A. R., 1974: Strong downslope windstorms. Mon. Wea. Rev. ,102, 592-602.

Black A., 1996: Personal communication. Department of Soil Sciences, University of British Columbia, Canada.

Durran, D. R., and J. B. Klemp, 1983: A compressible model for the simulation of moist mountain waves. Mon. Wea. Rev., 111, 2341-2361.

_____, D. R., 1986a: Another look at downslope windstorms, Part I, J. Atmos. Sci., 43, 2527-2543.

_____, D. R., 1986b: Another look at downslope windstorms, Part I, J. Atmos. Sci., 43, 2527-2543.

Eliassen, A., and E. Palm, 1960: Wave energy transfer in stationary gravity waves. Geofys. Publ., Vol. 22, Number 3.

Garratt, J. R., 1994: The atmospheric boundary layer. Cambridge University Press, 316pp.

Gutman, L. N., 1991: Downslope windstorms, Part I: Effect of air density decrease with height. J. Atmos. Sci., 48, 2545-2551.

Julian, L. T., and P. R. Julian, 1969: Boulders winds. Weatherwise, 22, 108-112.

Klemp, J. B., and D. K. Lilly, 1975: The dynamics of wave induced downslope winds. J. Atmos. Sci., 32, 320-339.

_____, J. B., and R. B. Wilhelmson, 1978: Numerical Simulation of convective storm dynamics. J. Atmos. Sci., 35, 1070-1096.

_____, J. B., and D. K. Lilly, 1979: The effects of terrain shape on nonlinearhydrostatic mountain waves. J. Fluid Mech., 95, 241-261.

_____, J. B., and D. R., Durran 1983: An upper radiation boundary condition permitting internal gravity wave radiation in mesoscale models. Mon. Wea. Review, Vol. 111, 430-444.

Lilly, D. K., 1983: Linear theory of internal gravity waves and mountain waves. Mesoscale Meteorology, D. Reidel Publishers, 781pp.

Long, R. R, 1953: Some aspects of the flow of stratified fluids, I: A theoretical investigation. Tellus, 5, 42-58.

Miller, P.P., and D. R. Durran, 1991: On the sensitivity of downslope windstorms to the asymmetry of the mountain profile. J. Atmos. Sci., 48, 1457-1473.

Raymond, D. J., 1972: Calculation of airflow over an arbitrary ridge including diabatic heating and cooling. J. Atmos. Sci., 29, 837-843.

Reisner, J. M., and P. K. Smolarkiewicz, 1994: Thermally forced low Froude numberflow past three-dimensional obstacles. J. Atmos. Sci., 51, 117-133.

Smith, R. B., 1980: Linear theory of stratified hydrostatic flow past an isolated mountain. Tellus, 32, 348-364.

Sullivan, P. P, J. C. McWilliams, and C. H. Moeng, 1994: A subgrid-scale model for large eddy simulations of planetary boundary layer flows. Boundary Layer Meteorology 71, 247-276.

Sun, W. Y., and C. Z. Chang, 1986: Diffusion model for a convective layer, Part 1: Numerical simulations of convective boundary layer. J. Climate and Appl. Meteo., 25, 1445-1453.

Whiteman, C. D., and J. G. Whiteman, 1974: An historical climatology of damaging downslope windstorms at Boulder, Colorado. NOAA Tech. Report ERL 336-APCL 35. [Available from the Superintendent of Documents, U.S. Government Printing Office, Washington D.C.]

Xue M. X, and S. J. Lin, 1991: On the numerical equivalence of advection terms in flux and advective form. Unpublished manuscript. Center for the Analysis and Prediction of Storms, University of Oklahoma, 10pp.