Since this is a new model, several tests were conducted to validate the numerical formulations, with five of the most important test groups presented here. The validation suite includes: two two-dimensional analytical mountain wave solutions, a simplified one-dimensional surface flux test using the Wangara data set, a three-dimensional density current simulation, and a two-dimensional simulation of the January 11, 1972 Boulder Colorado windstorm. With regard to the Boulder January 11, 1972 windstorm test, the current model is compared to established models. New results that are related to the proper numerical simulation of strong mountain waves are presented. Table 4.1 summarizes key model parameters used in each test.

Parameter	Linear Rigid	Linear Rad.	Long?s Equation	Wangara	Bubble	Boulder
nx	98	98	386	7	50	130
ny	4	4	4	7	50	4
nz	83	43	83	70	53	86
(m)	2000	2000	400	80000	800	1000
(m)	---	---	---	80000	800	---
(m)	200	200	250	30	100	341
(s)	20	20	10	20	4	5.0
(s)	5	5	1	20	1	2.5
(m/s)	20	20	10	0	0	Sounding
()	0.01956	0.01956	0.0108	Sounding	0	Sounding
()	250	250	---	Sounding	---	Sounding
()	--	---	---	---	-4.0	---
(m)	---	---	---	---	2000	---
(m)	---	---	---	---	14000	---
(m)	1	1	500*	0	500	2000
(m)	10000	10000	2000*	0	2000	10000
(m)	---	---	---	0	2000	---
	0.00001	0.00001	1.0e-5	0	5.0e-5	0.0004
	0.0	0.0	0.0	0	5.0e-5	0.00001
()	0.0025	---	---	---	---	0.0015
(m)	8000	---	---	---	---	18000
	0.2	0.2	0.2	0.2	0.2	0.2

*Terrain profile was numerically iterated using the non-linear lower boundary  condition, starting with an initial height of  = 570 m and width  = 2000 m.  Sounding data are listed in Appendix E.
 
 

4.1 Two Dimensional Linear Hydrostatic Mountain Wave

Two tests were performed using a linear hydrostatically forced mountain profile in an isothermal atmosphere. The goal here is to verify the model?s ability to reproduce basic mountain wave characteristics. In the first test, a Rayleigh sponge is applied to the top half of the domain (from 8-16 km) in combination with the rigid upper boundary condition (). The second test incorporates the upper linear hydrostatic  radiation condition at the 8km level. All other parameters remain unchanged and are given in Table 4.1.

The numerical solutions are compared to their analytical counterparts. Following Smith (1979), analytical solutions to a linear compressible hydrostatic mountain flow are developed in terms of streamline displacement from the undisturbed upstream value. For linearized compressible hydrostatic flow, the differential equation for streamline deflection is:

, (4.1)

.

The mountain profile is bell shaped:

                                                                      . (4.2)

,

where  is the air density at the surface () and  is the base state density. The analytical velocity fields are functions of the streamline deflection:

, (4.3)

. (4.4)

Results from the Rayleigh sponge rigid lid test are presented in Figures 4.1 and 4.2. Figure 4.1 displays the model, analytical, and difference perturbation horizontal velocities at a non-dimensional time  = 60 (30000 seconds). Figure 4.2 follows the convention of Figure 4.1 for the vertical velocity field. The present model reproduces the analytical fields with acceptable accuracy. The difference fields indicate a slight upstream tilt with height of the wave structure and accounts for most of the observed error. The maximum and minimum perturbations are within 10% of the analytical values at a height of one vertical wavelength (approximately 6.4 km). The upper radiation condition test requires only half the grid points in the vertical when compared to the sponge case. Figures 4.3 and 4.4 reveal similar characteristics found in the Rayleigh damping-rigid lid test. This suggests that for weakly forced hydrostatic flows the radiation condition can be used without appreciable loss in numerical accuracy. Both tests compare favorably with results presented by Nance and Durran (1994). They conducted a number of linear hydrostatic mountain wave experiments with similar initial and boundary conditions and computed difference fields to measure model accuracy. Their error characteristics are quite similar to the present models in terms of magnitude and location.

Another measure of the accuracy of a numerical model is the transport of wave energy with height. One type of wave activity, as discussed by Eliassen and Palm (1960), is the vertical flux of horizontal momentum. The momentum flux at the surface can be expressed by:

. (4.5)
 

(a)                                                                                                     (b)

 
 
 

                                                        (c)

 

                                                   (a)                                                                                      (b)

 

                                                (c)

 

Figure 4.2. Numerical model results for a linear hydrostatic mountain wave test using the Rayleigh damping rigid lid combination at time  = 60 (30000 seconds) for (a) , (b) , and (c) -. The contour interval for (a) and (b) is 0.0005 m/s, and (c) 0.000125 m/s. The area depicted follows that of Figure 4.1.
 

                                                        (a)                                                                                        (b)
  
 

                                                    (c)
 

Figure 4.3. Numerical model results for a linear hydrostatic mountain wave test using the upper radiation condition at time  = 60 (30000 seconds) for (a) , (b)  , (c) -. The contour interval is (a), (b) 0.005 m/s and (c) 0.001 m/s. The area depicted follows that given in Figure 4.1.

 
                                                      (a)                                                                                          (b)
  
 

                                                     (c)

Figure 4.4. Numerical model results for a linear hydrostatic mountain wave test using the upper radiation condition at time  = 60 (30000 seconds) for (a) , (b) , and (c) -. The contour interval in (a), (b) is 0.0005 m/s and (c) 0.000125 m/s. The area depicted follows that given in Figure 4.1.
 

The momentum flux at the surface is equal to but opposite in sign to the vertical flux of wave energy or surface wave drag, where and represent the perturbation velocities and  the base state density. The momentum flux is often normalized by the linear hydrostatic value.

,

where  is the surface reference density and  is the static stability for an isothermal atmosphere:

.

Figure 4.5 reveals a difference of only a few percent between the rigid lid/sponge and linear radiation condition tests. Both cases support momentum fluxes of 90% or better at one vertical wavelength and 95% of the normalized value at the surface. For steady state mountain waves, Eliassen and Palm (1960) show that the vertical flux of horizontal momentum remains constant with height when . The present model produces a nearly vertical profile and is consistent with other compressible models (Durran and Klemp, 1983 and ARPS, 1995). Normalized flux profiles approach unity for horizontal model domains twice that used above (not shown).
 

Figure 4.5. Vertical flux of horizontal momentum at  = 20, 40 ,60 normalized by the linear hydrostatic value. Solid lines represent the Rayleigh damping solution and the dashed lines represent the upper radiation condition results. One vertical wavelength is approximately 6.4 km.
 
 

4.2 Two-Dimensional Finite Amplitude Boussinesq Flow

The model is tested using a two-dimensional non-linear Boussinesq flow and compared to Long?s (1953) solution. Following Durran and Klemp?s (1983) implementation of Long?s solution, the differential equation for the streamline deflection from the undisturbed upstream height is similar to the linear case presented in Section 4.1. The equation for streamline deflection is:

. (4.6)

.

The Scorer parameter is . For a constant base state flow and static stability, Equation (4.6) describes non-linear flow using a linear differential equation. The lower boundary condition  is non-linear and cannot be directly applied. After implementation of the linearized upper radiation condition, the solution to (4.6) is

  .

A numerical approximation for the non-linear lower boundary streamline deflection is made using the initial terrain profile from (4.2). The procedure involves solving (4.6) for the surface streamline deflection using the initial terrain profile. The resulting streamline is substituted back into the non-linear lower boundary condition and iterated until the change in streamline deflection is within acceptable limits. Steffensen?s Algorithm found in Burden and Faires (1989) is used to obtain the iterated non-linear lower boundary condition for use in the numerical model. For an initial mountain height of 570 meters and mountain quarter-width of 2 km, the final mountain peak is 500 meters with the crest shifted upstream a few hundred meters. The analytically determined velocity perturbations have the same form as (4.3) and (4.4) without the base state density term:

, (4.7)

. (4.8)

The numerical model is prepared by applying the Boussinesq approximation to the buoyancy term in (3.3) and substituting a reference potential temperature into the total potential temperature for the pressure gradient terms of (3.1)-(3.3). The advection terms were removed from the pressure equation (3.5). In this form, the pressure is allowed to change according to the linearized divergence term. The input parameters for this test are provided in Table 4.1 under Long?s test group. The numerical solutions for are reproduced in Figure 4.6 at a non-dimensional time of = 60 (12000 seconds). Only the center portion of the domain (measuring 40 km across and 8 km deep) is shown. A comparison of the approximately steady state model predicted , , analytical , , and associated difference fields (- and-) is presented in Figure 4.7 and 4.8. Overall, the model reproduces the analytical solution reasonably well. The maxima and minima are within 7% for  and to within a few percent for  at a height of one vertical wavelength (approximately 5.8 km). The combined phase and amplitude errors for the velocities at the height of one vertical wavelength are on the order of 10-20%. Most of the difference is related to phase errors. Additional tests (not presented here) indicate that both the radiation and Rayleigh sponge/rigid lid applications produce similar errors in magnitude and phase. Note that these phase errors are consistent with those from the linear hydrostatic tests in Section 4.1.

For this test, the vertical flux of momentum is nearly constant with height and within a few percent of the normalized non-linear value at Ut/a = 60 (Figure 4.9). Analytical streamlines (trajectories involving a single time level of data) were computed and compared to the model equivalents with the results illustrated in Figure 4.10. The streamlines are computed by advancing a parcel through the respective velocity fields using a predictor-corrector method. The trajectory method and test solutions are presented in Appendix C. In general, the numerical model provides a good representation of the analytical streamline pattern over the mountain and displays errors characteristic similar to those shown in Figure 4.7 and 4.8.
 

                                                                            (a)

                                                                           (b)

  (a) (b)

Figure 4.6. Numerical solution for a non-linear Boussinesq forced mountain wave using the linear hydrostatic upper radiation condition for (a)and (b)  at time = 60 (12000 seconds). The contour intervals are (a) 0.25 K, and (b) 5.0 Pa. The plot is centered over the mountain with dimensions 40 km x 8 km.
 
 

                                                   (a)                                                                                              (b)

                                                   (c)

Figure 4.7. Numerical and analytical solutions and difference fields for the non-linear Boussinesq mountain wave using the linear hydrostatic upper radiation condition at time  = 60 (12000 seconds) for (a) , (b) , (c) -. The contour interval for (a) and (b) is 0.5 m/s and (c) 0.2 m/s. The plot is centered over the mountain with dimensions 40 km x 8 km.
 

                                                    (a)                                                                                            (b)

 

                                                (c)

 
Figure 4.8. Numerical and analytical solutions and difference fields for the non-linear Boussinesq mountain wave using the linear hydrostatic upper radiation condition at time  = 60 (12000 seconds) for (a), (b) , and (c) -. The contour interval for (a) and (b) is 0.25 m/s and (c) 0.1 m/s. The area shown follows that shown in Figure 4.7.

Figure 4.9. Vertical profile of the vertical flux of horizontal momentum for  = 15, 30, 45, and 60 normalized by the non-linear Boussinesq value. One vertical wavelength is approximately 5.8 km.
 
 

Figure 4.10. Comparison of numerically computed streamlines for the numerical model solution (solid line) and the analytical solution (dotted line) for the non-linear Boussinesq test (2) at a non-dimensional time  = 60.

4.2 Wangara Day 33

A simple test of the surface parameterization scheme is performed in conjunction with the Wangara Day 33 data set. The simulation is initialized with the observed 9 am vertical profiles of , , and . The primary focus of this test is to validate the sub-grid scale mixing of the potential temperature field when heat is applied to the surface. No horizontal gradients are present in the initial data set and therefore no gradients develop during the coarse of the simulation. Yamada and Mellor (1975) updated the horizontal wind fields with estimated geostrophic values computed from surface pressure data. Their results indicate that even with the estimated change in the geostrophic winds with time, the predicted profiles deviated from the observed data. Therefore, mixing of the horizontal wind components is not evaluated here. Since this is a one-dimensional test, it is an opportunity to tune the mixing length parameter for simulations not capable of resolving boundary layer convective cells. The relevant input parameters are listed in Table 4.1 under the Wangara test group. Following Sun and Chang (1986), heat flux was applied at the surface with a magnitude of 0.216 and half-period of 11 hours. Figure 4.11 compares the observed and model predicted potential temperature profiles as a function of time. The simulation compares favorably to the observed data. The turbulent mixing scheme captures both the average temperature and mixed layer depth.

                                                                                                Wangara Day 33 pt time/height

                                                                                                       Profile as Figure 4.11 (a)

                                                                                                                            (b)

                                               
Figure 4.11. Vertical  profiles from the (a) observed Wangara Test Day 33 from Sun and Chang, (1986), and (b) model at 9 am, 12 noon, 3 p.m., and 6 p.m. local time.
 
 

4.4 Three Dimensional Symmetry Test

ARPI3D is tested for symmetry errors by dropping a cold bubble over a symmetric mountain in a neutral environment. The initial bubble perturbation was -4 K with the mountain defined by (4.2). The model parameters used in this test are listed in Table 4.1 under the Bubble test group. This test is designed to show errors in the numerical implementation, which are asymmetric in nature. The horizontal velocities  and , are symmetric out to 11 decimal places at t=3000 seconds on a CRAY J90 computer. As the bubble moved out of the domain, slight asymmetries were noted (at the 10th place) in the horizontal velocity fields. This test was rerun using zero gradient lateral boundary conditions. The horizontal velocity fields remained symmetric out to the 14th place at t=3000 seconds. It appears that slight asymmetries are introduced by the lateral boundary condition scheme.

4.5 Idealized January 11, 1972 Boulder Colorado Windstorm Simulations

A number of researchers have numerically simulated the Boulder CO, January 11, 1972 windstorm event. The numerical results are often compared to the observed storm structures reported by Lilly and Zipser (1972) and Lilly (1978). Many of the observed windstorm characteristics, such as strong lee side surface winds and wave breaking regions in the upper troposphere and lower stratosphere, are well-represented by the numerical models. This is quite remarkable given the application of an idealized terrain profile (see equation 4.2), the estimated thermodynamic upstream profile, and the two-dimensionality of the numerical experiments.

This test group is broken up into three sections. The first section compares the present model?s numerical prediction to the historical numerical studies. The second part exposes the sensitivity of the numerical solution to the vertical resolution and advection scheme accuracy, and the third section addresses the lateral boundary condition issue.

4.5.1 Comparison with Previous Numerical Studies

Two well known two-dimensional numerical simulations of the Boulder January 11, 1972 windstorm were conducted by Peltier and Clark (1979) (denoted as PC) and Durran and Klemp (1983) (denoted as DK). Despite significant differences in terms of the system of equations, initial and boundary conditions, and solution techniques, results from the two models compare favorably. Both display amplification of the mountain wave and subsequent wave breaking structure in the lower stratosphere and upper troposphere with strong surface winds in the lee of the mountain crest. Numerical results from PC and DK and the present model, in terms of total potential temperature, are displayed in Figure 4.12 and 4.13. The present model?s results using a mountain shape defined by (4.2) and upstream conditions closely matches those reported by DK. Model parameters for this test are listed in Table 4.1 under the Boulder test group. The thermodynamic profile is evaluated from a sounding supplied by DK with the lowest pressure set to 820 mb and an isentropic layer between 820 mb and 685 mb. Data for levels above 110 mb were obtained from PC?s initial Grand Junction 1200Z sounding. For comparison with DK, the horizontal advection terms are approximated by a centered fourth order scheme. The vertical advection terms are computed using a centered second order scheme. The second order equivalent to that presented in (3.19)-(3.23) is obtained by setting the coefficients in the first terms to unity and omitting the second term. The observed cross section obtained from observations as reported by Lilly and Zipser (1972) is given in Figure 4.14. This model compares favorably with results from DK and PC. All three models generate similar transient gravity waves and wave breaking turbulent regions over the mountain crest.

As with previously reported aspects of the numerical storm, the wave response grows continuously and becomes unstable for the selected time step after approximately =10,000 seconds. A comparison of the surface wave drag as a function of time from DK, PC, and the current model are given in Figure 4.15 (a) and (b), respectively. The DK and PC models project an initial peak at approximately t=1000 seconds (due to the startup procedure) and then level off until about 6000 seconds, when the wave response and corresponding surface wave drag increases significantly. DK noted that the solution was unsteady and amplified with time, contrary to the observed storm, and suggests the unbounded growth may be due to the absence of surface friction. The solid line in Figure 4.15 (b) traces the surface wave drag computed from the present model. The surface drag time series from this model follows that given by DK and PC.

During the process of comparing the results from ARPI3D to previously reported numerical predictions, significant sensitivities were uncovered that were related to vertical grid resolution (or advection scheme accuracy), and lateral boundary conditions. These factors separate or combined, produce dramatic changes in the outcome of the experiments. The next two sections focus on these issues.
 

                                                        (a)                                                                            (b)                                                                (c)

Figure 4.12. Isentropes for the January 11, 1972 Boulder, CO windstorm at =4000s from (a) the present model, (b) DK, and (c) PC (at 4160s). The contour interval for (a) and (b) is 5K and unknown in (c).

 
                                                            (a)                                                                           (b)                                                                (c)

Figure 4.13. Isentropes for the January 11, 1972 Boulder, CO windstorm at =8000s for (a) the present model, (b) DK and (c) PC. The contour interval for (a) and (b) is 5K and is unknown in (c).

Figure 4.14. Observed potential temperature cross section from Lilly and Zipser (1972) for the January 11, 1972 Boulder Colorado windstorm.
 
 
 
 

                                                                                                                    (a)
 

                                                                                        Figure 4.15 (a) reproduced plot from

                                                                                            Durran and Klemp of surface drag...

                                                                                                                        (b)
                                 

Figure 4.15. Plots of surface wave drag as a function of time for (a) Durran and Klemp (solid line) and Peltier and Clark (dashed line), and (b) for ARPI3D.
 
 

4.5.2 Sensitivity to Vertical Resolution

In order to remove any potential lateral boundary effects and investigate the above noted sensitivities, two experiments were conducted in which the mountain is located at the 200-km mark of a 430-km wide domain. The experiments are distinguished by different vertical grid resolutions: dz=150m and dz=341m. All other model options and parameters remained unchanged except for a decrease in the time step to allow the simulation to continue during the amplification period without violating the CFL time step criteria. Both simulations, using second order vertical advection, place the wave breaking and turbulent zones at the base of the stratosphere, with the coarsely resolved case (dz=341m) exhibiting a stronger wave response in this region. A comparison of the total potential temperature fields is given in Figure 4.16. The wave activity in the stratosphere above the main breaking level is more coherent in the dz=150m case than in its coarser resolved counterpart. The vertical wavelength in the stratosphere is approximately 6km and resolved by 17 grid points in the coarse resolution and by nearly 40 grid points in the fine resolution case. A surface drag comparison for the vertical resolution tests is plotted in Figure 4.17. The coarsely resolved simulation displays a similar wave amplification pattern as the original simulation (Figure 4.15b), while the dz=150m case indicates an initial delay in the amplification process by approximately 5000s. The pseudo-steady state wave drag is noticeably lower for the dz=150m case. The amplification of the stratospheric momentum flux in the coarse resolution run is reduced (delayed), suggesting that the vertical transfer of wave energy and wave dissipation is sensitive to vertical grid resolution (Figure 4.18).
 
 

                                                                (a)                                                                                                                   (b)

Figure 4.16. Large domain (430x28km) model simulated isentropes for the January 11, 1972 Boulder, CO windstorm at =15,000 seconds for (a) dz=150m, and (b) dz=341m, after both cases achieved high drag states. Second order vertical advection is used in both runs. The contour interval is 10K and the entire model domain is plotted.

Figure 4.17. Surface wave drag for dz=150m (solid line) and dz=341m (dashed line) simulations. Note the dashed line follows the small domain results from Figure 4.15(b).

Figure 4.18. Vertical profile of the normalized vertical flux of horizontal momentum for the dz=150m case (solid line) and the dz=341m case (dashed line) at =5000 seconds. The flux was normalized by an estimated mid-tropospheric value. Note the increased wave activity near the tropopause for the dz=150m run.
 
 

4.5.3 Lateral Boundary Influence

As alluded to in the previous section, the lateral boundaries can pose a formidable threat to mountain wave prediction. In this section, two of the five boundary conditions discussed in Chapter 3 are tested for the dz=150m case. One method of testing lateral boundary conditions is first to perform a large domain control run in which boundary influences are minimal in a desired location. Then, conduct boundary condition sensitivity experiments for a smaller upstream domain, with the boundary in close proximity to the forcing mechanism. The control run for these tests is the high-resolution large domain run of the previous section. The two experiments involve the vertically averaged phase scheme of Durran and Klemp (1983) and the environmental phase speed. The environmental gravity wave phase speed was set to 2*. The choice of phase speed is somewhat arbitrary but it ensures that any gradient in the normal velocity component that reaches the boundary is passed on to the outside point. The inflow boundary (west boundary) is located 63.5km upstream of the mountain crest. This location corresponds to the 136.5km mark in the large domain control run. The simulations were advanced to =5,000 seconds and the perturbation  fields in the first 50km of the domain are compared in Figure 4.19 to the same area of the control run. In the control run, the perturbation fields are smooth, with significant wave energy present upstream of the mountain. Note that the velocity perturbations are nearly zero at the inflow boundary in the vertically averaged Orlanski case (Figure 4.19b). The environmentally determined phase speed compares rather well to the control run
 

                                                (a)                                                                                            (b)

 

                                                  (c)

 
Figure 4.19. Comparison of the perturbation horizontal velocity () at =5000 seconds for (a) large domain control run, (b) vertically averaged Orlanski phase speed, and (c) constant phase speed (). The contour interval is 2.0 m/s. The domain displayed is 50x28km and represents the same area with respect to the mountain peak in the control run (a). For (b) and (c) the left edge of the plot is the inflow boundary.
 

perturbation field, although there are still noticeable errors. The errors in the Orlanski case are of sufficient magnitude to hamper wave development over the mountain crest severely. In fact, the amplification stage present in the control run is delayed in the vertically averaged Orlanski run until approximately =20,000 seconds. The reason for this is tied to the development of the flow reversal region above the mountain. The inflow boundary reduces the negative  field in the vicinity of the tropopause to almost zero. The interior velocity field is unable to develop a wave-induced critical layer in a timely fashion. The amplification for the environmentally determined phase speed case occurs at approximately =9000 seconds, nearly 1 hour prior to the control run. The reason for this behavior is unclear, but may be due to the smaller domain size.

The formulation of the lateral boundary condition for the normal velocity component, following Orlanski, is inherently flawed. The advection equation approach in (3.17) and (3.18) requires the quantity () to be directed outward to modify the normal velocity at the outermost grid point. The computed phase speed is a function of the normal velocity field. As the solution approaches a steady state, the time tendency approaches zero, reducing the magnitude of the velocity used in the advection equation. For the vertically averaged Orlanski case, the phase speed and corresponding solutions for <2,000 seconds at the boundary is very similar to the control run. But as the solution advances, errors in the first order determination of the quantity () and spatial derivatives generate inward pointing phase velocities. The time dependent structure of the horizontal wind field at the boundary for the vertically averaged Orlanski case supports the findings of Durran et. al. (1993). Their study compares various phase speed estimation methods, including the Orlanski method, at the lateral boundaries in a single layer shallow water system. They reported that early in the simulation the Orlanski method reproduces the true phase speed reasonably well. But as partial reflections due to numerical errors occur at the boundary, the computed phase speed deteriorates significantly and becomes a poor representation of the true phase speed.
 

4.6 Model Test Summary

This newly developed numerical model compares favorably with a number of steady state analytical solutions and currently existing models in similar flow regimes. In regards to the well-known Boulder windstorm of January 11, 1972, the model performs admirably when measured against current models and observations. The windstorm tests reveal sensitivities not previously reported in terms of vertical resolution and lateral boundary condition type. Results from Sections 4.5.2 and 4.5.3 suggest that tighter constraints on the vertical grid resolution and lateral boundary conditions are needed in future windstorm prediction studies. Given the strong gradients of base state variables and the frequent development of a wave breaking region, the commonly used method of dividing the vertical wavelength by the grid spacing produces less than adequate results. Improved forecasts are possible by either increasing the vertical resolution or improving the accuracy of the advection scheme.