The Effects of Surface Friction on Downslope Wind and Mountain Waves

1 Department of Earth and Atmospheric Sciences, Purdue University, West Lafayette, IN 47907, USA 2 Department of Atmospheric Sciences, National Taiwan University, Taipei, Taiwan, ROC * Corresponding author address: Prof. Wu-Ron Hsu, Department of Atmospheric Sciences, National Taiwan University, Taipei, Taiwan, ROC; E-mail: hi@webmail.as.ntu.edu.tw A nonlinear nonhydrostatic model developed at National Taiwan University and Purdue University is applied to study the effects of surface friction on downslope wind and hydraulic jump for the windstorm of 11 January 1972 Boulder, Colorado, USA. The model reproduces strong downslope wind and hydraulic jump for both free-slip surface and viscous surface simulations. However, simulated hydraulic jump propagates downstream with an inviscid free-slip boundary, whilst the jump becomes stationary with a more realistic no-slip boundary condition. Furthermore, a returning flow (u < 0) and trapped waves also develop in the lower atmosphere downstream from the stationary hydraulic jump in the viscous surface case. These lee waves decay with distance due to surface friction. With the inclusion of surface friction, the model is able to generate a realistic turbulence field consistent qualitatively with aircraft observations at the time. (


INTRODUCTION
Mountain waves and downslope wind have been studied extensively since Queney (1947Queney ( , 1948)), Scorer (1949), and Long (1953aLong ( , b, 1954)), etc.Using the Grand Junction, Colorado sounding of 11 January 1972, many scientists (Klemp and Lilly 1975;Clark and Peltier 1977;Durran and Klemp 1983;Ikawa 1988;Doyle et al. 2000;Hsu and Sun 2001;Chen and Sun 2001) have applied numerical models to simulate the mountain waves, the hydraulic jump, and the strong downslope wind observed on the lee of the Rocky Mountains (Lilly and Zipser 1972;Lilly 1978).
Three different mechanisms have been proposed to account for the development of the strong downslope wind, including: (a) hydraulic theory-if the mountain height exceeds the threshold, a strong wind can develop along the lee when a subcritical flow transits to a supercritical flow (Long 1953a,b;Smith 1985;Durran 1986;Durran and Klemp 1987), (b) enhancement of lee-side surface wind by vertically propagating waves according to linear theory in an atmosphere with multilayer structure (Klemp and Lilly 1975;Lilly and Klemp 1979, etc.), and (c) downslope wind enhanced by the energy trapped by the wave-breaking region in the upper layer (Clark and Peltier 1977;Smith 1987).Detailed discussions are found in Durran (1990).
Most of the numerical studies on the subject are based on free-surface condition, which is partly due to the qualitative success of frictionless simulation in reproducing lee waves and downslope wind and partly due to the difficultly of introducing surface friction in a theoretical model (Richard et al. 1989).With the inclusion of surface friction, the numerical results obtained by Richard et al. (1989) eliminated the unrealistic high surface wind (reaching 70 ~ 80 m s 1 − or more) that propagates indefinitely downstream in free-slip simulations (Durran 1990).However, due to the very poor horizontal resolution used in their hydrostatic model, the hydraulic jump and the lee waves downstream from the jump are not well represented.Olafsson and Bougeault (1997) also showed diminishing wave breaking phenomenon in their no-slip simulation with a hydrostatic model.Hence, a fully nonhydrostatic NTU-Purdue Model (Hsu and Sun 2001) was used in this study.Our results not only confirm Richard et al's finding that the hydraulic jump becomes stationary on a viscous surface, but also reveal that wave trains develop in the lower atmosphere downstream from the hydraulic jump, which are not detected by their hydrostatic model.Other detailed mesoscale features, such as steep hydraulic jump, cloud formation, and distribution of the turbulence field are also discussed in the paper.

EQUATIONS
A simplified two-dimensional version of the NTU-Purdue Nonhydrostatic Model is used in this study.Our original model is three-dimensional and it is an extension of the Purdue hydrostatic Mesoscale Model (PMM) (Sun and Chern 1993;Chern 1994;Hsu and Sun 1994).The vertical coordinate of the NTU-Purdue Nonhydrostatic Model is defined as: where P 0 , pressure of a reference atmosphere, is strictly a function of height.Although this vertical coordinate appears to be a σ -pressure coordinate as used in many hydrostatic models, pressure in (1) is not a function of time.The position of each grid point is fixed in time.The model employs, strictly speaking, a σ -z coordinate.By ignoring the Coriolis force and y- component variation of all variables, the governing equations become: where Adv and Diff represent the advection and the diffusion operator, respectively.The advection terms are calculated with Sun (1993) forward scheme.The diffusion process is parameterized through a level 2.5 turbulence scheme.Density is a prognostic variable.The main advantage is that there is no diabatic term in the prognostic equation for density.Pressure field is diagnosed through equation of state, where T is temperature, and R is the gas constant.The equivalent potential temperature, θ e is used as the prognostic variable in the heat equation, where θ e is defined as: Here θ is potential temperature, c p is the specific heat at constant pressure, and q v is specific humidity of water vapor L v . is the latent heat of vaporization.With the definition of θ e , the effects of latent heat due to phase change do not appear explicitly in the first law of thermodynamic equation: The second term on the right hand side of equation ( 7) takes into account the change of θ /T.The prognostic equation for water in the model is: The total specific humidity, q q q w v l = + , is also a semi-conservative quantity in the absence of precipitation.Here q l , is the liquid water content.
We followed Deardorff's (1980) work for the turbulence scheme with the exception that all spatial derivatives are corrected for the terrain following coordinate.Eddy-coefficient relations were employed for fluxes of semi-conservative quantities, and the sub-grid scale eddy coefficient was made proportional to the sub-grid scale turbulence kinetic energy, E .The latter is governed by: where e u w ' ( ' ' ) ≡ + 1 2 2 2 and E e ≡ ' .The sub-grid fluxes were parameterized by: where K m is the sub-grid scale eddy coefficient for momentum, K h the sub-grid eddy coefficient for scalar quantities.The coefficients, A and B in (13) for unsaturated air are given by: and for saturated air, The sub-grid scale eddy coefficients in ( 13) ~ ( 17) were given by: where and l is a sub-grid scale mixing length which was required not to exceed the grid scale, ∆s or the height from surface, in magnitude.Following Sun and Ogura (1980), Deardorff (1980), the length scale in the stable atmosphere is assumed to be: Dissipation rate ε is The full turbulence parameterization is applied for both free-slip and no-slip simulations.As for the grid arrangement, the Arakawa C grid is employed.Density and w are also staggered in the vertical direction.As for time integration, the model uses a two-tier forward backward procedure, which is neutral in time with respect to both sound waves and internal gravity waves.In addition, the method is completely free from computation modes and very accurate.The detailed numerical scheme is presented in Hsu and Sun (2001).
As for the treatment of surface friction, the similarity equations (Businger et al. 1971) were applied between the ground surface (where u = v = w = 0) and the layer at z = 25 m for the simulation with surface friction.It is also assumed that neither potential temperature gradient nor heat flux exchange exists between soil and the atmosphere.When frictional effects are omitted, a free-slip boundary condition is applied at the lower boundary (i.e., ∂ ∂

NUMERICAL RESULTS
The domain consists of 471 grids in the x-direction with ∆ x = 1 km and 400 levels in the z-direction with ∆ z = 75 m, except the first level above the ground surface where ∆ z = 25 m.The mountain profile is a witch of Agnesi curve, z x h x a s ( ) /( / ) = + 1 2 2 , where h =2 km and a = 10 km, which are typical values for the Colorado Front Range (Doyle et al. 2000).The initial condition (same for both free-slip and no-slip boundaries) is horizontally homogeneous and based upon the upstream 1200UTC 11 January 1972 Grand Junction, Colorado, sounding shown in Fig. 1.The same initial condition has been used in Doyle et al. (2000).
(a) Free-slip case with h = 2 km.
The virtual potential temperature, x-component wind u, vertical velocity w, and perturbation pressure (departed from horizontal average) after 2.5 h of integration within the domain of 140 km ≤ x ≤ 300 km and 0 ≤ z ≤ 18 km are shown in Fig. 2. The simulated downslope wind and vertical velocity reach 70 m s 1 − and 25 m s 1 − , respectively.They are comparable with those presented in Hsu and Sun (2001), Chen and Sun (2001), and others discussed in Doyle et al. (2000).It is noted that a semi-implicit scheme is used to filter acoustic waves in Chen and Sun (2001), which is quite different from the model presented here.The enhancement of u-component wind in Fig. 2b is consistent with the pressure perturbation from the horizontal average shown in Fig. 2d.A local maximum for both wind and pressure perturbation occurs at x = 220 km near surface, and another local maximum at z ≅ 9 km and x = 240 km, which are also consistent with Chen (1999).
The simulation after 3.5 h and the time sequence of vertical velocity at z = 3625 m after 1.5, 2.5, and 3.5 h are presented in Fig. 3. Results show that the downslope wind reaches 80 m s 1 − at the surface; both hydraulic jump and strong downslope wind propagate farther downstream with time.In addition to gravity waves and over-turnings occurring in the upper atmosphere, a few branches of vertical motion develop in the lower atmosphere between the mountain peak and hydraulic jump that propagates farther downstream and allows more space for the development of organized vertical motions, as shown in Figs.3c, f.Figures 3a, d also   show that cloud forms approximately between two constant virtual potential temperature (θ v ) lines of 295K and 310K, in both upstream and downstream regions of the hydraulic jump.
There is no cloud between the mountain peak and hydraulic jump, where the air is relatively warm and dry, even though the vertical motion is strong.
Figure 3e shows the distribution of Richardson number.It is calculated through the negative ratio of the buoyancy production and shear production terms in equation ( 9) of the previous section.This ratio (flux Richardson number) is then multiplied by K K m h / to obtain the gradient Richardson number shown in Fig. 3e.The shaded area (with positive buoyancy) in the middle troposphere above the lee-side slope of the mountain corresponds to the overturning of airflow induced by the strong downslope wind.The air is also very turbulent downstream from the hydraulic jump (x > 250 km) in the lower troposphere.It is noted that the propagating hydraulic jump also acts as an obstacle to the westerly flow in the upper level just as the mountain does to the airflow in the lower atmosphere.The jump induces wave breaking in the layer approximately 10 km < z < 15 km, and to the east of the jump, x > 250 km.The air to the upstream of the jump (i.e., x < 250 km) at this level, on the other hand, is free of turbulence.

(b) Viscous surface case with h = 2 km
With the inclusion of the surface frictional effect, the model generates a downslope wind of about 50 m s 1 − near the ground surface and a hydraulic jump with vertical velocity of 25 m s 1 − after 2.5 h of integration, as shown in u, v, θ v , and perturbation pressure (Fig. 4).The maximum velocities are weaker than the inviscid case as in Fig. 2 because of the surface friction and a shorter distance between the mountain peak and the hydraulic jump for the downslope flow.Note that wave trains with significant amplitude show up downstream from the hydraulic jump.This feature is seen in the aircraft observations at the time (Fig. 7 of Lilly 1978).The wave amplitude and wavelength decrease with the increase of distance.Furthermore, the returning flow (u < 0) forms in the boundary layer to the east of the jump.This feature can be seen more clearly in more detailed structures of the downslope wind and hydraulic jump presented in Fig. 4e.Peng and Thompson (2003) used the distribution of the Richardson number to identify the boundary layer depth in their study of the surface frictional effect on flows over mountains.The distribution of the Richardson number is also plotted in Fig. 4e and it revealed that the boundary depth is very shallow for airflow over a high mountain, as the air becomes turbulent (small Richardson number) very close to the ground surface, especially on the lee side slope of the mountain.This result is consistent with Peng and Thompson (2003).
Comparing the simulated results after 2.5 h (Fig. 4) and 3.5 h (Fig. 5), we see that the location of hydraulic jump does not move with time on a viscous surface, which is consistent with Richard et al. (1989).Consequently, the distribution of clouds (Fig. 5d) and Richardson number, R i (Fig. 5e), are different from the inviscid case.Meanwhile, the horizontal length scale of the waves to the east of hydraulic jump increases with time (Figs.5f and g).The train of lee waves is apparent in the potential temperature and vertical velocity at z = 3625 m after 1.5 h (dotted line), 2.5 h (dashed line), and 3.5 h (solid line).The amplitude of those waves is about 5 m s 1 − or more in the vertical motion; and a few degrees in virtual potential temperature.Eventually the flow pattern in the lower atmosphere becomes quite different from the case with a free-slip surface condition.Figure 5h show the Reynolds stress at the surface.The surface Reynolds stress can be expressed as the square of the friction velocity (u*).We can examine the effect of surface friction through the following momentum equation: Since the jump remains in the same location (nearly steady) and the vertical motion is comparatively weak near the surface, the du/dt term can be approximated by In addition, the diffusion term near the surface is dominated by the surface frictional effect.Thus, where δh is the height of the surface layer.The magnitudes of each term in the momentum equation near the surface can then be estimated by: The mathematical operator, ∆ , takes the horizontal difference of a variable over a distance L, which is the distance for the downslope wind braking to a full stop near the ground surface.Since the maximum intensity of the surface wind is 50 m s 1 − at x = 210 km, and u is 0 at x = 215 km (Fig. 4b), du/dt is approximately -0.25 m s −2 .The pressure gradient force term is approximately -0.14 m s −2 , with ρ ~1 kgm −3 , ∆p~700 Pa (Fig. 4d), L~5 km.The last term of the equation is about -0.1 m s −2 , with u * 2 ~0.2 m s 2 −2 (Fig. 5h), δh~20 m.These three terms roughly come to a balance and the analysis shows that the surface friction plays a key role in stopping the downslope wind.Surface friction consumes the energy of the downslope wind so that the location of hydraulic jump remains the same.The stationary hydraulic jump also acts as an obstacle to the westerly flow to produce gravity waves to the lee side, which is similar to mountain waves triggered by dryline or surface front reported by Sun and Wu (1992).In the inviscid case, the areas of strong disturbances spread far away from the mountain, because of the propagation of hydraulic jump.As for the simulation with surface friction, strong disturbances are confined closer to the mountain downstream from the jump, which is more consistent with the observation (Compare Fig. 5e and Fig. 7 of Lilly 1978).The stationary jump also triggers upper level waves that are locked approximately in the same location.Thus, the phase remains the same and the upper level waves (Fig. 5c) become much stronger than those (Fig. 3c) generated by the propagating jump in the simulation with a free-slip ground surface.
It should be noted that although our result is consistent with Richard et al. (1989) in terms of the stationary property of the jump for the frictional case, the hydraulic jump produced by Richard et al. (1989) must be very weak with the use of a very poor horizontal resolution (10 km) in their hydrostatic model (the magnitude of the upward motion not reported in the paper).In addition, the mountain waves are stacked vertically right over the top of the mountain ridge, instead of tilting toward the downstream direction.Furthermore, the lee waves downstream from the jump are completely missing in their simulation due to the poor spatial resolution.As for our simulation, we are able to generate a jump much more realistically with a nonhydrostatic model.The jump is very strong with a maximum vertical velocity of 25 m s 1 − and it still maintains a stationary property as in Richard et al. (1989) with the inclusion of surface friction.
The mountain waves propagate in the downstream direction and the simulated lee waves are consistent with the observation presented by Lilly (1978).However, the downslope wind in our results is weaker than that obtained by Richard et al. (1989) due to the differences in the initial wind profiles in the two studies.Richard et al. (1989) used the initial maximum wind speed of 50 m s 1 − at z~10 km, proposed by Peltier and Clark (1979).We used an observed   3 except for the simulation with a viscous surface after 3.5 h of integration.Panel (g) is the virtual potential temperature, θ v , at z = 3625 m and (h) surface Renolds stress ( u * 2 ) after 1.5 h (dotted line), 2.5 h (dashed line), and 3.5 h (solid line).profile upstream from the Boulder Windstorm with a maximum wind of 40 m s 1 − occurred at z~7 km (Fig. 1).Hence, they obtained a stronger downslope wind (60 m s 1 − ) than ours (50 m s 1 − ) for the simulation including surface friction.

(c) Simulations for h = 500 m
As the mountain height decreases, the intensity of downslope wind and hydraulic jump decrease for both ground surface conditions.The wave trains become the dominant feature for mountain height of 1 km or less (Durran 1986).Our results also show that hydraulic jump disappears when the mountain height is 500 m (Fig. 6), although waves and turbulent mixing still exist in the middle and upper atmosphere in our simulations.It is noted that near the mountain peak, the upward motion reaches only 0.9 m s 1 − (Fig. 6a) for the simulation on a free-slip surface.The simulation with surface friction produces a stronger upward motion with maximum over 1.2 m s 1 − (Fig. 6b).This phenomenon can be explained by examining Fig. 7.
Since the wind speed reduces to zero at the slope of the mountain for the simulation with surface friction, the magnitude of the horizontal divergence in the lower atmosphere (z < 2.5 km) is larger than that of the free-slip simulation.However, the waves dissipate more quickly over a viscous surface than a free-slip surface.

SUMMARY
The NTU-Purdue Nonhydrostatic Model was applied to study the effects of surface friction on downslope wind and hydraulic jump for the Boulder windstorm of 11 January 1972.With a free-slip surface, the model generates a strong downdraft wind reaching 80 m s 1 − and upward motion of 25 m s 1 − after 3.5 h of integration.Meanwhile, the downslope wind, hydraulic jump, and associated waves in the upper atmosphere propagate downstream with time.
On the other hand, the model with a viscous surface produces weaker and more realistic downslope wind.Surface friction consumes the energy of the downslope wind so that the location of hydraulic jump remains the same.The hydraulic jump becomes stationary, which is consistent with Richard et al. (1989).In addition, wave trains develop downstream from the hydraulic jump, which is also consistent with aircraft observations shown by Lilly (1978).The model also generated trapped wave trains for mountain height of 500 m.Because of surface friction, the excited waves were stronger in the beginning, but decayed faster on a viscous ground surface.

Fig. 1 .
Fig. 1.Vertical profiles of the (a) virtual potential temperature (K) and (b) u-component ( m s 1 − ) from the Grand Junction, Colorado, USA sounding for 1200UTC 11 January 1972.

Fig. 3 .
Fig. 3. Simulated (a) virtual potential temperature, θ v (with a contour interval of 10K), (b) u (10 m s 1 − ), (c) w (5 m s 1 − ), (d) liquid water content, g l (0.0001), and (e) Richardson number, R i after 3.5 h; panel (f) the vertical velocity, w, at z = 3625 m after 1.5 h (dotted line), 2.5 h (dashed line), and 3.5 h (solid line) with a free-slip surface and 2 km mountain.The areas with a negative Richardson number are shaded in panel (e); and only contour lines of R i = 0.3 and 0 are drawn.

Fig. 4 .
Fig. 4. Panels (a)-(d) are the same as Fig. 2 except for the simulation with a viscous surface after 2.5 h of integration.Panel (e) is the zoom-in view of the airflow near the vicinity of the mountain.The magnitude of the vertical motion is amplified by 5 times to show the jump motion.The distribution of Richardson number is also plotted.The areas with a negative Richardson number are shaded; and only contour lines of R i = 0.3 and 0 are drawn.

Fig. 5 .
Fig. 5. (a)-(f) are the same as Fig.3except for the simulation with a viscous surface after 3.5 h of integration.Panel (g) is the virtual potential temperature, θ v , at z = 3625 m and (h) surface Renolds stress ( u * 2 ) after 1.5 h (dotted line), 2.5 h (dashed line), and 3.5 h (solid line).

Fig. 6 .
Fig. 6.Simulated vertical velocity, w, after 3.5 h of integration (a) with a freeslip surface and (b) with a viscous surface.The maximum height of the mountain is 500 m for these two experiments.Contour interval is 0.3 m s 1 −

Fig. 7 .
Fig. 7. Simulated horizontal velocity, u, and corresponding horizontal divergence field, ∂ ∂ u x , after 3.5 h of integration (a, b) with a free-slip surface and (c, d) with a viscous surface.The maximum height of the mountain is 500 m for these two experiments.Contour interval is 2 m s 1 − for u, and 0.0002 s −1 for ∂ ∂ u x .