GRAVITY WAVE DRAG PARAMETERIZATION
1. Introduction

Studies of systematic errors show NWP models have a tendency toward a positive bias in the equator-to-pole pressure gradient and too strong a westerly wind in the midlatitude northern hemisphere (Boer, et.al., 1985, Palmer, et.al., 1985: PSS). The importance of sub-grid scale momentum transfer to the global momentum budget of the atmosphere from the breakdown of mountain induced gravity waves was realized by Bretherton (1969), Lilly (1972), Lindzen (1981) and others. Observational studies of Lilly, et.al., 1982) and Brown (1983) and theoretical modeling studies of Eliassen and Palm (1960), Peltier and Clark (1979, 1983) suggest that the momentum flux divergence from sub-grid scales can affect the motions of the large scale circulation. PSS and Miller and Palmer (1987) found that if the residual of the zonal and time-averaged momentum budget is calculated from analyses, then the magnitude of the residual term of the zonal-mean momentum budget is significant with largest negative contribution over latitudes where orography is prominent. This indicates the possibility of a missing sub-grid scale drag term in the model's momentum equation. Enhancing orography, spatially varying the roughness and/or drag coefficient as well as directly compensating for forecast errors are momentum damping mechanisms that have been tested to correct model deficiencies described above, but each of them suffer from drawbacks. While these effects usually show up on climatological time scales of GCM experiments, NWP models have improved their medium range performance skill by utilizing gravity wave drag as a momentum damping mechanism.

The gravity wave drag parameterization used at NMC is a fairly crude approximation of gravity waves forced by flow over topography. The surface flux parameterization follows Pierrehumbert (1987). With respect to the vertical deposition of momentum, we favor the saturation flux method (Lindzen, 1981) of depositing the momentum vertically after studying comparisons with a linear deposition. As described below, critical layers, the degree of instability and other processes contribute to the vertical momentum deposition within the saturation flux method. An attempt is made to include these contributions.

Basic Formulation

The scheme described here is taken from the ideas of Pierrehumbert (1986) (PH) and those of Helfand et. al. (1987) who followed the work of PSS and Lindzen (1981). The surface stress is calculated using linear theory for small Froude number. For large Froude number the surface stress calculation attempts to incorporate non-linear effects which act to limit the growth of the momentum flux. The philosophy used is that a momentum flux source is at or near the ground and has a strength determined by the variance in the earth's elevation, low level wind and atmospheric stability. This gives rise to a "base level" momentum flux, which is the net amount that may be applied to the layers above. The stress is distributed conservatively in the vertical as a function of a number of criteria.

In interior layers of the atmosphere, if no absorption occurs, the stress is constant with height and no momentum deposition occurs. At the highest sigma layer, any remaining stress is deposited. Otherwise the stress is not absorbed with height until either 1) a critical level is encountered at which point the momentum flux vanishes or 2) the gravity waves break and generate turbulence from shear or convective instability or 3) from wave saturation. Once the gravity wave begins to break the wave saturation hypothesis of Lindzen (1981) is invoked. That is, the waveinduced turbulent dissipation prevents the displacement amplitude from exceeding its critical value, in this case 1 - 1/(4Ri), where Ri is the Richardson number. The Froude number, Fr, and the Richardson number of the mean state of the atmosphere, Ri, are related to the wave modified Richardson number, Rm, as Rm = Ri - Ri Fr. If the critical value of the wave-modified Richardson number is chosen to be 1/4, the above criteria results. This is discussed in more detail below, in terms of application to the NMC spectral model.

The drag due to sub-grid scale momentum flux appears as an additional body force on the right hand side of the momentum equation of the model. The starting point is the momentum eqn from page . The velocity tendency equation is modified by the vertical derivative of the gravity wave induced stress:

where is surface pressure.

The prognostic equations for momentum, in the NMC spectral model, (see Sela, 1982, 1987) are in the form of tendencies of vorticity and divergence. The relevant part for this section can be represented as:

where A and B are the velocity tendencies and include the vertical diffusion of momentum. The vertical flux of momentum due to gravity wave drag is incorporated in the model by modifying A and B to become:

in eqn. (6.2).

The vertical derivative of the stress is written in difference form as

where SI(k) is the index for model sigma interface levels and SL(k) is the index for the model sigma layers.

Based on linear, two-dimensional non-rotating, stably stratified flow over a ridge, gravity wave motions are set up which propagate away from the mountain. The flux measured from a low level in the atmosphere is defined to be the base level flux. It is taken to be the first 1/3 of the atmosphere but this choice is arbitrary. The vertical momentum flux or gravity wave stress in a grid box due to a single mountain is given as in PH:

where is a grid increment, N is the Brunt Viasala frequency given as

The environmental variables are calculated from a mass weighted vertical average over the base layer. In eqn. (6.5) G(Fr) is a monotonically increasing function of Froude number, , where U is the wind speed calculated as a mass weighted vertical average in the base layer, and , is the vertical displacement caused by the orography variance. An effective mountain length, where m is the number of mountains in a grid box, can then be defined to obtain the form of the base level stress used in the computation as shown in eqn. (6.7).

giving the stress in a grid box. PH gives the form for the function G as

where is an order unity non-dimensional saturation flux set to 1.0 and 'a' is a function of the mountain aspect ratio. Typical values of U=10m/s, N=, = 100km, and a = 1 give a flux of 1 Pascal. If this flux goes to zero linearly with height, then the deceleration would be about 10m/s, which is consistent with numerous observations (see Table 1 in PH). In Fig. 6.1, the universal flux function, G, is shown as a function of Froude number and mountain shape parameters . A value of 1.0 for 'a' has been used in the experiments following.

When breaking occurs for some layer in the free atmosphere, the stress at the next lower layer is known, beginning with the base level flux described in (6.7). It should be noted that this representation of the stress (6.7), under a binomial expansion for small Froude number, is identical to that used in the free atmosphere in PSS:

with , where is an inverse length scale parameter given as in the PSS formulation.

Following Helfand et. al. (1987) and PSS wave breaking occurs at the level where the amplification of the wave causes the local Froude number to exceed the critical value

where is the critical vertical displacement amplitude for a wave in the free atmosphere which is about to break. This means that a shear instability is assumed to occur at Ri<1/4. PH points out that a criteria for shear instability is difficult to determine and convective breaking (Ri<0) criterion is usually satisfied before that of shear breaking. Both of the above criteria are used in the present parameterization. However, if Ri<0 occurs in the base level, then no deposition takes place.

Once it is determined that some wave breaking is to occur, the vertical structure of the stress is reduced so that the wave-induced turbulent dissipation prevents the displacement amplitude, , from exceeding its critical value. This is measured by the ratio of the critical and locally calculated Froude numbers described above. The gravity wave momentum flux should be reduced by a factor when the Froude number exceeds its critical value. Using (6.9) we can rewrite the Froude number at interface in terms of the stress of the interface layer below the one under consideration, , as

As a result the vertical profile of the momentum flux is built up from the base level stress as a function of local values of N, U and the particle displacement. The gravity wave drag-induced velocity tendency change is computed as the vertical divergence of the momentum flux as described above in eqns. (6.3, 6.4).

Program Implementation and Coding Considerations

There is an option to treat the first KBPS-1 sigma levels as the interpolation between sigma level KBPS and the base layer (bottom stress) when calculating the vertical stress profile. For example, if KBPS is set to 2, then sigma layer 1 is not calculated by (6.11) but by a linear interpolation between layer 2 and the base layer stress. This is done so that the first sigma layer values are governed by the model planetary boundary layer formulation and the thickness of the first sigma level is less then the mountain variance.

Another option independently allows for a linear interpolation of the vertical stress profile to begin at an arbitrary layer LCAP and extend to the model top sigma layer. This option allows a linear profile to be imposed in part or entirely throughout the atmosphere. Linear profiles of the vertical stress function have been used with success on the NMC Nested Grid Model, and at other NWP centers.

The scheme is composed of subroutines, V3GWD, where basic values and base level stress are calculated, as well as the base level stress and velocity tendencies, and GWDRAG which receives the input from V3GWD and returns the vertical structure of the stress. The changes to the velocity tendency are calculated from the vertical stress profile calculated in GWDRAG. Within the base layer, no gravity wave deposition is allowed if a critical layer or instability takes place. The base layer is assigned as the layers below 667mb. The parameter KBPS indicates the sigma level where the saturation flux hypothesis is to be used. Below this level a linear stress profile is fit to the bottom stress value. The current default value for KBPS is 2. The parameter LCAP is used to replace the saturation flux hypothesis with a linear profile extending from sigma layer LCAP to the top of the model atmosphere. Changes to the control variables NUM and CON may be used to alter the parameterization (see Chapter 8). Fig. 6.2 shows sigma interfaces and levels of the current 18 layer model. The base level averages of variables on sigma layers and interfaces have been treated separately. The result is that stress is defined on an interface while the velocity tendency is in sigma layers. The variables, symbols in the subroutine calls are shown in Table 6.1.

Table 6.1

List of variables and symbols passed between the subroutines.

Sigma coordinate approximate pressure (mb) Interface Layer
SI(k) SL(k)

19 _______________ 000

(KDIM) 18 --------------- T, u, v, A, B 021

18 _______________ , Ri, Vg, Fr 050

. . . .

KSM _______________ , Ri, Vg, Fr

KSM --------------- T, u, v, A, B

. . . .

4 _______________ , Ri, Vg, Fr 948

3 --------------- T, u, v, A, B 960

3 _______________ , Ri, Vg, Fr 973

2 --------------- T, u, v, A, B 981 (KBPS)

2 _______________ , Ri, Vg, Fr 990

1 --------------- T, u, v, A, B 995

1 _______________ , Ri, Vg, Fr 1000

Figure 6.2: Model vertical structure as it applies to the gravity wave drag parameterization

References

Brown, P. R. A., 1983: Aircraft measurements of mountain waves and their associated momentum flux over the British Isles, Quart. J. Roy. Meteor. Soc 109, 849-865.

Boer, G. J., N. A. McFarlane, R. Laprise, K. D. Henderson and J. P.Blanchet, 1984: The Canadian Climate Centre Spectral Atmospheric General Circulation Model, Atmos. Ocean., 22, 397-429.

Bretherton, F. P., 1969: Momentum transport by gravity waves, Quart. J. Roy. Meteor.Soc, 95, 213-243.

Brown, P.R.A., 1983: Aircraft measurements of mountain waves and their associated momentum flux over the British Isles, Quart. J. Roy. Meteor. Soc , 109, 849-865.

Eliassen, A., and E. Palm, 1960: On the transfer of of energy in stationary mountain waves, Geophys. Publ., 17, 1-44.

Helfand, H.M., J.C. Jusem, J. Pfaendtner, J. Tenenbaum and E. Kalnay, 1987: The Effect of a gravity wave drag parameterization scheme on GLA fourth order GCM forecasts. Submitted to JAS.

Lilly, D.K., 1972: Wave momentum flux: A GARP problem, Bull. Amer. Meteor. Soc., 20, 17.

________, J. M. Nicholls, R. M. Chervin, P. J. Kennedy and J. B. Klemp,1982: Aircraft measurements of wave momentum flux over the Colorado Rocky Mountains, Quart. J. Roy. Met. Soc., 108, 625-642.

________, Lindzen, R. S. 1981: Turbulence and stress due to gravity wave and tidal breakdown, J. Geophys. Res., 86, 9707-9714.

Miller, M.J. and T.N. Palmer, 1987: Orographic gravity-wave drag: its parameterization and influence in general circulation and numerical weather prediction models. Submitted to Quart. J. Royal Met. Soc.

Palmer, T.N., G.J. Shutts and R. Swinbank, 1986: Alleviation of a systematic westerly bias in general circulation and numerical weather prediction models through an orographic gravity wave drag parameterization, Quart. J. Roy. Meteor. Soc., 112, 1001-1039.

Peltier, W.R., and T.L. Clark, 1979: The evolution and stability of finite-amplitude mountain waves. Part II: Surface wave drag and severe downslope windstorms, J. Atmos. Sci., 36, 1498-1529.

_________, and ________, 1983: Nonlinear mountain waves in two and three spatial dimensions, Quart. J. Roy. Meteor. Soc., 109, 527-548.

Pierrehumbert, R.T., 1987: An essay on the parameterization of orographic gravity wave drag. Geophysical Fluid Dynamics Laboratory/N.O.A.A., Princeton University, Princeton, NJ 08542

Sela, J., 1982: The NMC spectral model. NOAA Technical report NWS 30, U.S. Department of Commerce, N.O.A.A., National Weather Service, 36pp.

________, 1987: The new T80 NMC operational spectral model. Preprint, Eighth Conference on Numerical Weather Prediction, Baltimore, Maryland, February 22-26, 1988.