A new three dimensional model (with a two dimensional option) is used to simulate dry heated mountain wave flows. The model, ARPI3D, was constructed initially in a simple two-dimensional framework to test horizontal boundary condition applications for use in the ARPS. The code is extensively documented and the initialization of the horizontally homogeneous base state variables follows that of the ARPS.

Since its inception, ARPI3D was modified to include a simple coordinate transformation, a vertically-implicit solving technique, and a linearized upper radiation condition between pressure and vertical velocity. In addition, this model solves for the non-dimensional Exner function  instead of pressure. Initially, ARPI3D included a "p" (pressure) or "pi" (Exner function) option but the p-option was later dropped for the more efficient pi-system. The inclusion of the non-linear pressure gradient and divergence term coefficients closes the model energy budget. The decision to include the total potential temperature in the pressure gradient term is not solely based on the goal of energy conservation. It was due in part from tests comparing the linearized and total term versions. Perturbations on the order of 30% of the base state potential temperature were observed in the breaking wave regions in the stratosphere.

3.1 Equation Set

The ARPI3D equation set is:

  , (3.1)

  , (3.2)

  , (3.3)

  , (3.4)

  . (3.5)

The total Cartesian velocities , , and  are defined as:

 ,  is the base state wind in the x direction.

 ,  is the base state wind in the y direction.

 , is the velocity in the z direction.

The thermodynamic variables, non-dimensional pressure, potential temperature, and their base state equivalents are computed according to:

 , and , with  ,

.

The coriolis terms  and  are:

,

,

where  is the angular velocity of the Earth and  the latitude. Equations (3.1)-(3.5) describe the Euler equations with the addition of a heat source  in the potential temperature and pressure equations and turbulent and computational mixing terms in (3.1)-(3.4). Normally, mesoscale models make a number of approximations that are focused 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 and are incorporated in the current model configuration. The non-linear divergence terms is included to close the model energetics. The non-linear pressure gradient term reduces the time step slightly but improves the numerical solution. Little improvement is noted when the coefficient of the divergence term is at full strength. Equations (3.1)-(3.5) conserve energy in the absence of friction, mixing, and heat.

The above equation set can be transformed to a terrain following coordinate system by:

,

where is the top of the model domain, is the height of the surface, and  is the physical height of the computational surface. Using the chain rule, the coordinate transformation is given by:

,  ,  .

The spatial derivatives in the transformed system become:

  , (3.6)

  , (3.7)
 
  , (3.8)

 

  , (3.9)

  , (3.10)

  , (3.11)

  , (3.12)

  . (3.13)

The contravariant vertical velocity is defined by,

,

or following (3.6)-(3.8),

.

The base state variables  are horizontally homogeneous and hydrostatically balanced.

.

The turbulent mixing terms are described in detail in Section 3.2. The numerical smoothing terms  are composed of computational mixing and Rayleigh damping:

,

,

,

.

Sections 3.4.1 and 3.4.2 give a detailed description of the computational mixing and Rayleigh damping terms.
 
 

3.2 Sub-Grid Scale Closure

The sub-grid scale processes are those that cannot be predicted explicitly by the model. In the planetary boundary layer (PBL), the sub-grid scale mixing can contribute significantly to the evolving flow field. Thus, considerable work has focused on estimating the unresolved flow patterns. When the surface is heated and the atmosphere becomes unstable, convection occurs initially on scales proportional to the mixed layer height. There are a number of ways in which to represent convective processes in the boundary layer by the model. Two types of mixing mechanisms are presented and encoded in the model.

The first method applies strong diffusion in areas of neutral or unstable lapse rates. This sub-grid closure scheme requires large eddy viscosities capable of mixing the unstable boundary layer and developing the height of the mixed layer in a timely manner. Since the convective elements are not explicitly resolved, this type of scheme allows for large horizontal and vertical grid spacings on the order of 5000m and 200m, respectively. This method is preferred if the convective motions in the boundary layer are small when compared to the surrounding flow.

A second method of addressing the unstable boundary layer is the Large Eddy Simulation (or LES). The purpose of the LES is to explicitly resolve the majority of the convection in the heated boundary layer. As shown by Deardorff (1980) and many others, this requires resolutions on the order of 100 meters in each spatial direction. This approach makes use of the sub-grid closure scheme, but the amount of energy in the unresolved scales is significantly reduced. Most researchers are satisfied when 90-95% of the fluxes, energy, and variances are explicitly resolved.

For large three-dimensional mountain wave flows, the LES resolution requirements place a severe constraint on the size of the domain that can be studied with current computing power. Horizontal resolution on the order of 400-1000 m is suitable for simulating mountain waves but is not sufficient to resolve properly the convective elements in the mixed layer. Thus, simulations of broad mountain ranges are better served by the eddy viscosity method than the LES approach. The scale analysis of Chapter 2 indicates that strong mountain wave flows are more energetic than boundary layer circulations. But in the interest of evaluating the sub-grid closures, both methods are tested and the results are compared in Chapter 5.

The turbulent parameterization scheme used in this study follows the base model of Sullivan et. al. (1994). It was first developed by Deardorff (1980) and later applied by Moeng (1984). This scheme is chosen because of its simplicity and ability to match the simulated planetary boundary layer solution to Monin-Obukov similarity theory in a LES. The method utilizes a 1.5 order closure scheme that predicts turbulent kinetic energy () for use in the determination of the eddy viscosity. For simulations with horizontal resolutions greater than 200m, the convective motions in the early development stage of the boundary layer cannot be adequately resolved. The mixing length in the convective boundary layer is enhanced following Sun and Chang (1986). This modification was incorporated after tests the using the Wangara Day 33 data set and the baseline Deardorff eddy viscosity model predicted an untimely development of the mixed layer height. The unmodified Deardorff and Sun and Chang schemes are presented later in this section. The sub-grid closure terms are given in Cartesian coordinates and equated to their coordinate transformed counterparts via:

,

,

,

.

The strain tensors are defined by:

,

,

,

,

,

.

By definition , , . Following Deardorff (1980) the eddy viscosity is:

,

where  = 0.1 . The velocity stresses () in Cartesian and transformed coordinates are given by:

,

,

,

,

,

.

The scalar fluxes are defined by:

,

,

.

The eddy diffusion is defined in terms of the mixing length, average grid spacing, and the eddy viscosity using:

.

The average grid spacing is defined by the relation . The sub-grid scale kinetic energy (), used in defining the mixing length  and the eddy viscosities  and , is predicted following:

,

where  is the shear production term,  the buoyant production term,  the diffusion of  , and  the dissipation of . The sub-grid scale energy source or sinks terms are defined as:

,

,

,

,

where  = 0.93 and each index  are summed from 1 to 3. The mixing length above the boundary layer is computed following Deardorff (1980):

for neutral and unstable stratification ,

, for stable stratification ,

where  is the grid computed static stability.

In the Sun and Chang application, Deardorff?s length scale in the unstable boundary layer is replaced by a fraction of the observed peak wavelength for the vertical velocity component (). Caughey and Palmer (1979) obtained an expression for the peak  for experimental data collected in the Kansas, Minnesota, and Ashchurch field studies. The following expression provides a good fit to the observed data:

,

,

where is the depth of the mixed layer and z the height above the ground surface. The last term inside the square bracket is small and is here omitted. The inversion height is diagnosed by comparing the surface potential temperature to the potential temperature in and above the mixed layer.
 
 

3.3 Surface Flux Parameterization

A simple surface drag formulation is used to represent the surface stress. The method presented here follows that implemented by Miller and Durran (1991) and can be found in Haltiner and Williams (1980). The surface stress for the  and  velocity components is computed by:

, (3.14)

. (3.15)

,

,

where  and  are the base state surface wind components. This formulation reduces the stress for only positive perturbations. The negative perturbations remain unchanged. The drag coefficient is computed from the relation:

,

where is the height at which the drag coefficient is valid. The above is either substituted in for the flux in the explicit method or is a source term in the implicit method. The surface heat flux is introduced as a flux or a volume source term. The source term has the form:

  , (3.16)

3.4 Numerical Smoothing

3.4.1 Computational Mixing

Two types of numerical smoothing are used in this study. The first type, computational, acts to remove small-scale structures created by non-linear aliasing and dispersive effects created by the advection scheme. This smoother is of fourth order and is applied to  in the horizontal and vertical directions. The form follows that used in the ARPS and is computed in computational space using:

,

,

,

.

The specified mixing coefficients are  and , and are generally chosen to be on the order of times the grid spacing to the fourth power.

3.4.2 Rayleigh Damping

Additional smoothing is used to damp gravity waves in the upper part of the modeling domain when the top boundary is set to the rigid lid condition. This type of numerical damping is designed to simulate a radiation condition by preventing reflection of gravity wave energy off of the top boundary (rigid lid). For this type of damping to be effective, the damping layer needs to be greater that one vertical wavelength in depth and a minimum of 30 grid points (Klemp and Lilly, 1978 and Durran and Klemp, 1983). Commonly referred to as Rayleigh damping, it is applied to  and is defined here by:

, for .

The damping terms are:

,

,

,

,

where is the height at which the damping begins and  is the height of the top of the model. The coefficient  represents the maximum damping coefficient () at the top of the damping layer, generally selected so that the dominant horizontal wave number is damped after about 25 time steps.
 
 

3.5 Boundary Conditions

Boundary conditions can pose a formidable challenge to long term mesoscale numerical predictions and to simulations that interact significantly with flow adjacent to the boundaries. The only true physical boundary in this numerical model is the bottom boundary. Much effort has been spent developing lateral and vertical boundary conditions for models whose domains are not periodic. To my knowledge, a fully robust open lateral boundary condition, which properly handles a combination of acoustic, gravity, and inertial waves, does not exist. As a result, this model includes a variety of schemes for predicting the time dependent variables at the lateral and vertical boundaries.
 
 
 

3.5.1 Lateral Formulation

There is a need to specify the normal advection and velocity terms at the first grid point outside the physical boundary. The normal velocity boundary condition is presented first. ARPI3D includes five choices for computing the boundary values of the normal velocity at the lateral boundaries. The reason for such a variety is simple: each specific modeling application could require a different type of boundary condition. The five lateral boundary condition options for the normal velocity components are:

a) the zero gradient condition

b) the Orlanski (1976) condition

c) the vertically averaged Orlanski phase speed scheme of Klemp and Lilly (1979), later modified by Durran and Klemp (1983)

d) the Klemp-Wilhelmson (1978) constant phase speed method

e) a newly developed hydrostatically (or environmentally) estimated phase speed
 

Schemes (b)-(e) are designed to allow waves in the interior of the model domain to pass freely out through the horizontal boundary with minimal reflection. With the exception of the zero gradient condition (which is self-explanatory), all of the above methods estimate the gravity wave phase speed and replace the normal velocity component equation of motion with outflow advection:

  , (3.17)

 

  . (3.18)

 

For method (d), the gravity wave phase speeds  and  are specified by estimating the fastest gravity wave phase speed. For method (e), the phase speed is estimated by the linear hydrostatic value. It can be found from the relation:

For hydrostatic modes  and we have:

                     ,

where the vertical wave number  is the hydrostatic limit corresponding to the fastest wave. The method is also applied to the  velocity component. For moist convection, Klemp and Lilly (1979) found it was advantageous to over-estimate the phase speed rather than underestimate it.

A review of (c), from the Durran and Klemp (1983) implementation, is given here. The vertically averaged phase speeds are formulated by solving (3.17) and (3.18) for the local gravity wave phase speed and averaging only the outgoing components. The inward directed components are set to zero prior to averaging. The vertically averaged phase speeds  for  and  are computed using the following relationships:

  ,

  ,
 
  ,

  ,

For methods (b) and (c) the phase speed is computed following Orlanski (1976). Note that the original Orlanski form is obtained by removing the vertical averaging step in the Durran and Klemp procedure.

For the horizontal advection terms at the lateral boundary, the higher order advection schemes cannot be directly applied. At an outflow boundary, the stable first order one-sided upstream differencing scheme is used for the normal horizontal advection term. On inflow boundaries, the normal advective term is set to zero since there is no information on the gradient of the advected variable outside the boundary.

The lateral boundary conditions for the fourth order computational mixing terms follows those used in the ARPS. The intermediate computational mixing value (second order derivative) at the point just inside the boundaries is used to set the value of the point outside the boundary. All other quantities that require definition of intermediate quantities at the boundary computations, such as the stress and strain terms in the sub-grid scale closure model, invoke the zero gradient condition.
 

3.5.1 Vertical Formulation

The commonly used Rayleigh sponge-rigid upper boundary and fixed bottom boundary conditions are employed in this model. In addition, following Klemp and Durran (1983), a radiation condition between  and  is available at the top boundary. The radiation condition forces the numerical solution at k=nz-1 for  and k=nz-2 for  with a linear analytical Boussinesq relation. The analytical solution allows upward propagating gravity wave energy to pass vertically out of the top boundary and removes the need for a sponge in the top part of the model domain. This method can reduce significantly the vertical extent of the computational domain without a degradation of the solution in a wide variety of gravity wave problems. The radiation condition is given later in this section.

When the rigid boundary condition is applied to the top and bottom boundaries, the perturbation scalar quantities in the vertical, with the exception of  at the surface, are set to zero gradient following:

 

, and  ,

, and  ,

, and  .

.

The vertical velocity and the contravariant vertical velocity at the lower boundary are defined by the relationship between the slope and the horizontal wind speed and the impermeability condition, respectively. The vertical velocity at the lower boundary is

.

The contravariant vertical velocity is

.

For the rigid upper boundary condition, both vertical velocities are set to zero

.

Quantities involved in intermediate steps at the boundary for the turbulent and computational mixing terms make use of the zero gradient condition, regardless of the top boundary condition type.

The upper radiation condition makes use of the linear Boussinesq hydrostatic system of equations. This formulation follows that developed by Klemp and Durran (1983):

,

where and  are Fourier transformed non-dimensional pressure and vertical velocity,  is the local static stability, and are the horizontal wave numbers,  is the local base state potential temperature, and is the specific heat. A comparison of the upper radiation and Rayleigh sponge-rigid lid upper boundary conditions for a simulated linear hydrostatic mountain wave is given in Chapter 4.

3.6 Discretized Equations

The model equations are computed using spatially centered finite differences on the Arakawa (1966) staggered C-grid. The grid box is illustrated in Figure 3.1. The velocity variables are located on the sides of a grid box and define the physical model boundaries. The scalar quantities are defined at the center of each grid box. The model scalar variables are computed on scalar points from 1 to nx-1 and 1 to ny-1 in the horizontal plane and from 2 to nz-2 in the vertical direction (Figure 3.2). The  velocity is computed from 1 to nx , 1 to ny-1, and 2 to nz-2 . The  velocity component

Figure 3.1. ARPI3D grid box displaying the spatial arrangement of the scalar and vector quantities.

Figure 3.2. Location of the scalar and velocity points in terms of each of the model axes. Hatching indicates the physical model domain.

is computed from 1 to nx-1, 1 to ny, and 2 to nz-2. The vertical velocity is computed from 1 to nx-1, 1 to ny-1, and 2 to nz-1.

The spatial derivatives for the gravity, inertial, and acoustic wave terms are differenced using centered second order accurate finite differences. The time integration follows the split time step approach of Klemp and Wilhelmson (1978). The time integration is broken up into two parts: terms evaluated using big time step information, denoted by the superscript  corresponding to advective, gravity, and inertial wave modes as well as all mixing terms and sub-grid scale processes. The terms evaluated on the small time step, denoted by the superscript , represent the acoustic modes and include the pressure gradient and divergence terms. The forcing terms computed on the big time step are used in the small time step to advance the dependent variables forward in time. This approach has been implemented and tested by a number of modeling efforts including Durran and Klemp (1983) and the ARPS (1995). The full discretized equation set is:

,

,

                                                                                     ,

,

                                                                             .

Note that all the variables except for  are advanced on the small time step using the forward time scheme. The forcing terms for  and are applied in the small time step. The small time step is advanced using the forward scheme in  steps where

                                                                                                                  .

In the non-dimensional pressure equation, the scheme is referred to as forward-backward differencing since the updated velocities are used in the divergence term. The potential temperature is evaluated and integrated on the big time step using leapfrog time differencing. The non-dimensional pressure and  are advanced using the fully explicit forward scheme or the vertically implicit forward scheme. The implicit scheme removes the vertical grid spacing time step restriction for acoustic wave propagation. The vertically implicit scheme follows the Crank-Nicholson method and is absolutely stable with respect to vertical sound wave propagation. The Crank-Nicholson scheme requires more computations due to the need to solve a tridiagonal matrix, but is beneficial for horizontal to vertical grid spacing ratios greater than 2.5. Use of the 2 time step mode splitting method reviewed above allows for certain unstable acoustic modes to exist. Durran and Klemp (1983) suggest that values of  between 0.5 and 1.0 successfully damps the unstable modes. These unstable modes are due to interactions between the sound waves and the advective and gravity waves as shown by Skamarock and Klemp (1992). A description of the vertically implicit solving technique is given in Appendix A. Appendix B reviews  upper boundary condition used in the vertically implicit time marching scheme.

The spatial averaging and differencing operators are defined by:

,

.

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). This form conserves the first and second order moments.

                                         ,  (3.19)
                                         ,   (3.20)
                                        ,  (3.21)

                                 , (3.22)

The contravariant vertical velocity is computed using

.

The sub-grid scale contributions are discretized according to:

,

,

,

.

The stresses are computed using:

,

,

,

,

,

.

The sub-grid scale kinetic energy is integrated on the big time step using the leapfrog scheme.

The advection and shear production terms are:

,

.

The strains are:

,

,

,

,

,

.

The mixing terms  are represented by:

,

,

,

.

The computational mixing coefficients  and  are held constant and generally chosen to be on the order of 0.0005. Finally, an Asselin (1972) is applied to all the predicted variables on the big time step to prevent the divergence of the odd-even time step solutions:

.