| Home > Newsevents > Training > Rcourse_notes > DATA_ASSIMILATION > GRAV-WAVE_CONTROL > |
The control of gravity waves
in data assimilation
27 April 1999
By Adrian
Simmons
European Centre for Medium-Range Weather
Forecasts
Table of contents >>
Next Section >>
Previous Section >>
2 . Non-linear normal-mode initialization
2.1 Basic method
Consider a forecast model linearized about a state of rest which is statically stable. The governing equations of the model can be written in the following general form:
|
(1) |
where
is a column vector of suitably-scaled model
variables, of dimension 
say, and 
is a symmetric, real 
matrix. Explicit
examples of sub-elements of 
will be given in the following subsections. The eigenvalues of 
are the frequencies of the small-amplitude wave
motions that the resting state can support, and the corresponding eigenfunctions
describe the structures of the waves. In limiting cases these waves can
be separated into a set of "meteorological" or "Rossby" waves which move
relatively slowly westward, and sets of eastward- and westward-moving gravity
modes, the graver of which move much faster then the Rossby waves. More
generally, the categorization is complicated in the tropics by the existence
of the so-called mixed Rossby-gravity and Kelvin waves. Comprehensive studies
of the modes of the shallow-water equations (known as the "Hough modes")
have been made by Longuet-Higgins(1968) and Kasahara(1976).The general non-linear model equations can be written in the form:
|
(2) |
where the term
represents both the non-linear terms and the residual linear terms
that arise from differences between the actual state of rest and that used
to compute 
.Consider a subset comprising
orthonormal
eigenvectors, 
, of 
with large eigenvalues 
. These represent the gravity waves that are to be initialized. A
general model state 
can be written in the form:
|
(3) |
Under idealized conditions in which the
eigenvectors represent a complete set of gravity waves, 
represents the component of the atmospheric state comprising
solely Rossby waves. More generally, 
represents the residual state orthogonal to the set of 
gravity modes, which comprises Rossby waves, the slower gravity waves
that will not be initialized, and mixed waves. We denote the component that
will be initialized by 
:
|
(4) |
The
are determined by taking the scalar products
of the two sides of (3) with the 
, and using orthogonality, 
, and orthonormality, 
, to obtain
|
(5) |
or
|
(6) |
where
is the 
matrix whose columns are the 
, for 
.Substituting (3) into (2) and projecting onto the
modes to be initialized gives:
|
(7) |
where
.2.1 (a) Linear normal-mode initialization:
Linear normal-mode initialization comprises simply setting to zero the amplitudes
of the 
modes to be initialized. It thus involves replacing
the uninitialized state 
by 
, where
|
(8) |
and
is the identity matrix.Fig. 7 presents an example of the impact of linear normal-mode initialization on the evolution of surface pressure at two points, one in the Great Plains of North America and one in the Himalayas. It is taken from the first study of normal-mode initialization for the original ECMWF finite-difference model (Temperton and Williamson, 1981; Williamson and Temperton, 1981), and presents an extreme case in that the initial data had been taken from an analysis produced using a different model. Imbalance in the interpolated initial conditions for the ECMWF model resulted in large-amplitude oscillations in the surface pressure, as indicated by the solid lines in the figure. The amplitude of these oscillations was clearly reduced by applying linear normal-mode initialization (dashed lines), but significant fluctuations remained. Also worthy of note is the large, rapid fall in surface pressure at the Himalayan point at the start of the integration from the initialized analysis.
Figure 7 . Surface pressure (hPa) as a function of time before (solid) and after (dashed) linear normal-mode initialization at 40oN 90oW (left) and 30oN 90oE (right), from Temperton and Williamson(1981).
2.1 (b) Non-linear normal-mode initialization:
Linear normal-mode does not prevent gravity waves from growing to significant amplitude because the non-linear or residual linear terms represented by
are not negligible. In non-linear normal-mode initialization
(Machenhauer, 1977), the initialized gravity-wave
amplitudes 
are set to non-zero values which are chosen to make the gravity-wave
tendencies small. An approximate solution of the equation
|
(9) |
is sought. Equation (7) is used to define an iterative solution of (9):
 .
|
(10) |
Two iterations of the procedure are usually found to be sufficient. The starting values could be taken from a linear normal-mode initialization (
), but in practice it is found to be
sufficient to start from the uninitialized analysis (
). This is the most straightforward starting point when the procedure
is implemented in terms of updates applied directly to model variables,
as outlined below in equations (12) to (15). It is, moreover, the necessary starting
point for an effective full-field diabatic initialization using the approach
described in 2.5 (a).An example of the working of non-linear normal-mode initialization is given in Fig. 8 . This is for the same case as presented in Fig. 7 for linear normal-mode initialization. Fig. 8 also contains information related to the vertical decomposition of normal modes which is discussed in the following subsection. For the moment, attention should be directed first towards the upper two panels and the curves labelled "0" which show the forecast from uninitialized initial conditions and "5" which show the forecast from the conditions most fully initialized by the non-linear normal-mode technique. These two panels can be compared directly with Fig. 7 which shows the forecast from linear normal-mode initialization. It is clear that the non-linear technique is much more successful than the linear technique. It not only prevents high-frequency oscillations in surface pressure over the 24-hour integration period, but also avoids the large change to the initial surface pressure at the Himalayan point, a change that is rejected by the model in the first few time steps of the forecast from linear initialization. Comparison of the lower-right and upper-left panels of Fig. 8 shows how the first iteration is sufficient to remove most of the short-period oscillations present when running from the uninitialized analysis, with the second iteration removing most of what is left after the first iteration.
Figure 8 . Surface pressure (hPa) as a function of time at 40oN 90oW (upper left) and 30oN 90oE (upper right) with no initialization (solid) and after non-linear normal-mode initialization of the first three and five vertical modes with two iterations. Surface pressure is also shown at 40oN 90oW with one, two and four modes initialized with two iterations (lower left), and with five modes initialized with one and two iterations (lower right). From Williamson and Temperton(1981).
Equation (10) is not normally solved directly. (7) and (10) are combined to express the iterative non-linear normal-mode initialization procedure in terms of the modifications made to the gravity-wave amplitudes at each iteration:
|
(11) |
Equation (11) is then converted into the initialization increment applied to model variables. Starting the iteration from the uninitialized analysis, the non-linear initialization procedure becomes:
|
(12) |
|
(13) |
where
is the 
diagonal matrix of eigenvalues 
for 
. The non-linear model tendency 
is computed using the latest approximation, 
, to the initialized state:
|
(14) |
The final initialized state,
, is given by:
|
(15) |
where
is the number of iterations.This "explicit" form of non-linear normal-mode initialization thus involves a repeated set of operations comprising:
| • One timestep of the model to computes the tendencies in physical space, starting from the latest approximation to the linearized state; |
| • A projection of these tendencies onto the gravity modes to be initialized; |
| • Division of the gravity-mode tendencies by the frequencies of the modes; |
| • Projection back to physical space; |
| • Addition of the increment to form a new approximation to the initialized state. |
Equation (13) may be written in the alternative form:
|
(16) |
Equation (16) is used as the starting point for the derivation of "implicit" non-linear normal-mode initialization. As discussed in subsection 2.4, this approach uses a modified form of the governing, linearized equations which enables increments to be calculated in physical space without use of explicit projections to and from gravity-wave space.
2.2 Vertical decomposition
Normal-mode initialization is a practical method because it is possible to separate the vertical and horizontal dependence of the gravity-wave modes. Furthermore, it has to be applied only to a limited number of modes in the vertical.
We consider small-amplitude perturbations about a state of rest with temperature
and surface pressure 
, with eastward and northward wind perturbations 
and 
, and temperature perturbation 
. The equations are set out in the form appropriate for a model (such
as that of ECMWF) in which the logarithm of surface pressure is a basic
prognostic variable, with perturbation given by 
. We restrict the presentation here to the case in which the
reference temperature 
is isothermal, rather than varying with pressure. This simplifies
the equations, and is not a serious restriction in practice. Effects of
unrepresentivity of the reference state of rest are taken into account in
non-linear normal-mode initialization through the term 
. Moreover, gravity-wave phase speeds
and vertical mode structures are not strongly dependent on the choice of
reference state, as illustrated in Appendix A.The dry, linearized, primitive equations for a terrain-following vertical coordinate are then:
|
(17) |
|
(18) |
|
(19) |
|
(20) |
Here
is the divergence:
|
(21) |
is the pressure-coordinate vertical velocity:
|
(22) |
and
is the temperature-dependent part of the perturbation
geopotential:
|
(23) |
is the planetary rotation rate, 
is the planetary radius, 
is longitude, 
is latitude and 
is the gas constant of dry air. 
, where 
is the specific heat of dry air at constant
pressure. In evaluating the pressure on a coordinate surface in equations
(19), (20),
(22) and (23), the surface pressure is taken to be
the reference value 
.We consider a model with
levels, define 
to be the column vector comprising the
values of 
at the 
levels, and define 
, 
and 
similarly. The column vector with each element equal
to 
is denoted by 
, and we use matrices 
and 
, and vector 
to represent a general finite-difference scheme:
|
|
|
|
Explicit expressions for
, 
and 
are given in Appendix A
for the vertical discretization scheme developed for the ECMWF model by
Simmons and Burridge(1981).We further define:
|
(24) |
Then, using (17), (18) and the combination of (19) and (20), the linearized equations become:
|
(25) |
|
(26) |
|
(27) |
where
|
(28) |
The normal-mode initialization procedure determines an increment to the "mass" variable
. To obtain
increments in temperature and the logarithm of surface pressure, equation
(27) is used together with
the vector forms of the equations for 
and 
:
|
(29) |
|
(30) |
Given initialization increment
, the corresponding increments to temperature and the logarithm of
surface pressure, 
and 
, are given by
|
(31) |
and
|
(32) |
Now let
denote the diagonal matrix of eigenvalues
of 
, for 
, and let 
denote the matrix whose columns are the eigenvectors of 
. 
can be written: 
.If the prognostic variables are transformed in the following way:
|
|
and
|
equations (25), (26) and (27) become a set of
uncoupled shallow water equations with "equivalent" depths 
:
|
(33) |
|
(34) |
|
(35) |
Here
is the divergence associated with velocities
and 
as defined
by (21):
|
(36) |
The governing equation for pure gravity waves is obtained by setting the rotation rate,
, to zero,
taking the time derivative of (35), and using (33) and (34) to eliminate the rate of change of divergence
. This gives:
|
(37) |
For scales small enough for a local plane geometry to be valid, the phase speed of the pure gravity waves is
. On the sphere the pure gravity-wave modes are the spherical harmonic
functions 
, where 
is an associated Legendre function. The spherical harmonic functions
are eigenfunctions of the Laplacian operator with eigenvalues 
. The pure gravity-wave frequencies are thus 
.Examples of modes are presented here for the current operational 31-level ECMWF model and for a 50-level version of the model planned for operational implementation. The location of the model levels for these two vertical resolutions is shown in Fig. 9 . Table 1 shows phase speeds and equivalent depths of the faster modes for a reference temperature,
, of 300K and a reference surface pressure, 
, of 1000hPa. Corresponding vertical structures of the divergence
fields (weighted by 
) are shown in Fig. 10 .
Figure 9 . The distribution of full model levels for 31-level (left) and 50-level (right) vertical resolutions, plotted for a distribution of surface pressure which varies from 1013.25 to 500 hPa.
Mode number  |
Phase speed (ms-1)for 31-level model |
Phase speed (ms-1)for 50-level model |
Equivalent depth (km)for 31-level model |
Equivalent depth (km)for 50-level model |
|---|---|---|---|---|
| 1 |
343 |
347 |
12.0 |
12.3 |
| 2 |
203 |
262 |
4.2 |
7.0 |
| 3 |
119 |
194 |
1.4 |
3.8 |
| 4 |
78 |
145 |
0.6 |
2.1 |
| 5 |
56 |
112 |
0.3 |
1.3 |
| 6 |
90 |
0.8 |
||
| 7 |
75 |
0.6 |
||
| 8 |
63 |
0.4 |
||
| 9 |
54 |
0.3 |
Figure 10 . Structures of the vertical modes of the 31-level and 50-level models with gravity-wave phase speeds faster than 50 ms-1, for a reference temperature of 300K and a reference surface pressure of 1000 hPa.
The gravest mode moves the fastest, and differs in structure from the other modes in that it is the only mode whose
-weighted amplitude (or energy density) decreases with height. It
is generally referred to as the external gravity wave in the context of
numerical weather prediction, but is known in geophysical fluid dynamics
as the Lamb wave. It is shown in Appendix B that there is a corresponding
analytical solution of the continuous equations in which the divergence
and temperature vary in the vertical as 
, and the phase speed is given by 
. This gives a phase speed of 347ms-1 for a reference
temperature of 300K, with 
=287m2K-1s-2
and 
=0.286. This value is reproduced to the nearest ms-1 by
the 50-level resolution.The remaining ("internal") modes provide a mathematically-acceptable basis for representing the vertical structure of model variables, but do not correspond to modes of the continuous equations. They have an oscillatory structure with amplitude approximately proportional to
, and represent standing waves that can exist because of the reflective
nature of the upper boundary condition 
applied in the model (Lindzen et al., 1968).It can be seen from Table 1 that only a limited number of modes are associated with gravity-wave phase speeds of the order of 50ms-1 or more. It is found necessary to initialize only these modes. Indeed, applying the non-linear normal-mode technique to the slower, higher-order modes is counter-productive, as the iterative solution fails to converge for these modes. The example presented in Fig. 8 illustrates the extent to which high frequency oscillations are suppressed by initializing up to the five gravest modes. Initializing four or five modes appears adequate in this case. These results were for a nine-level vertical resolution for which the gravity-wave phase speeds of the fourth and fifth modes were 70ms-1 and 39ms-1 respectively, for a 300K reference temperature.
2.3 Horizontal decomposition
The horizontal decomposition is illustrated most conveniently for the case of a global spectral model such as used at ECMWF. Firstly, the shallow-water momentum equations are rewritten in terms of the divergence,
(given by (36)), and
relative vorticity, 
:
|
(38) |
The non-linear shallow-water equations are:
|
(39) |
|
(40) |
and
|
(41) |
These equations are recast in terms of the spectral coefficients
, 
and 
in spherical-harmonic expansions:
|
(42) |
|
(43) |
and
|
(44) |
As the prognostic variables are real valued,
, 
and 
are the complex conjugates of 
and 
, and we need consider only the equations for 
. These equations can be separated into 
sets, one for each zonal wavenumber 
. We introduce scaled variables:
|
(45) |
|
(46) |
and
|
(47) |
For each
, we use 
, 
and 
to denote column vectors of dimension 
with elements 
, 
and 
, for 
.The shallow water equations (39), (40) and (41) may then be written in the form:
|
(48) |
where
and 
are diagonal matrices, with diagonal elements 
and 
, where
|
(49) |
and
|
(50) |
is a symmetric tridiagonal matrix with diagonal
elements zero, and off-diagonal elements 
for 
, with
|
(51) |
In the absence of rotation
, and the non-trivial eigenvectors of 
are gravity waves with frequencies 
, as noted earlier.With rotation, the basic Coriolis effect (as occurs in an f-plane geometry) is represented by the matrix
, while the matrix 
represents the variation of the Coriolis
effect with latitude, the so-called "beta" effect. Pure barotropic Rossby
waves are obtained by imposing 
and have frequencies 
.Equation (48) provides an explicit form for
:
|
(52) |
Eigenvalues and eigenfunctions of
can readily be obtained using standard software for matrix analysis.
Separation into Rossby and gravity modes can in practice be achieved simply
on the basis of wave frequencies; the eastward moving waves and faster westward
moving waves are the modes to be initialized.2.4 Implicit initialization
Although formally straightforward for global spectral models, explicit normal-mode initialization has a demanding requirement for computer storage of the normal modes in the case of high horizontal resolution. Similar considerations apply to global models based on alternative horizontal discretizations. Two additional problems may arise for a limited-area model. Firstly, the map projection may cause horizontal separability to be lost. Secondly, there may be a problem in defining lateral boundary conditions for the normal modes. This led to the development of implicit normal mode initialization. This is presented here for the case of a global spectral model; Temperton(1988) introduced the method for a regional model using a polar stereographic projection and finite elements with a non-uniform grid.
We start from the formula for initialization increments given by equation (16):
|
The approach relies on being able to determine the "fast" gravity-wave component
from the total tendency 
and on being able to solve the linear system for the increments 
. This can be achieved by approximating the operator 
. The approximation is to neglect the "beta-effect" in the
vorticity equation. Equation (52) is replaced by:
|
(53) |
Equation (49) indicates that the approximation will be a good one for small horizontal scales (large
) for which the non-zero elements 
of the diagonal matrix 
are small.With this approximation, the linearized form of (48) gives:
|
(54) |
|
(55) |
and
|
(56) |
The Rossby modes are stationary and non-divergent:
|
(57) |
and, from (55), satisfy the balance equation:
|
(58) |
From (54) and (56), a gravity mode (which has
) must have:
|
(59) |
The expanded form of (16) is
|
(60) |
Temperton(1989) showed that (57), (58) and (59) can be used to derive formulae for initialization increments that do not require explicit computation of gravity modes:
|
(61) |
and
|
(62) |
with
|
(63) |
The matrix
is pentadiagonal, but the equations
to be solved ((61) and (62)) can in fact each be split into two
separate equations for odd and even spectral components with diagonally-dominant
tridiagonal matrices on the left-hand sides. Their solution is thus straightforward.In practice this implicit method has been found to give results that are negligibly different from explicit normal-mode initialization for most of the wavenumber range of high resolution models. A mixed scheme has been adopted for the ECMWF spectral model. The implicit method is used by default to initialize total wavenumbers in the range
, with the explicit method used for 
, as the neglect of the 
term becomes increasingly inaccurate as the total wavenumber decreases.2.5 Diabatic initialization
The formalism of non-linear normal-mode initialization presented in 2.1 is general, but in practice it was found necessary to suppress the diabatic (parametrized) components of the term
in early applications of the method (Williamson and Temperton, 1981). In particular,
the large and rapidly varying temperature tendencies arising from the parametrization
of convection were found to inhibit convergence of the method. Moreover,
the basic initialization condition expressed by (9):
|
is inappropriate for the westward-moving thermal tidal waves forced by the daily westward progression of the solar heating of the atmosphere.
The form of initialization in which only the adiabatic tendencies from the model are computed at each iteration is referred to as adiabatic non-linear normal-mode initialization. Although successful in preventing the growth of high-frequency oscillations in subsequent forecasts, adiabatic initialization also suppresses the slowly-varying large-scale circulations that are a feature of the tropical atmosphere. Two of the approaches that have been adopted to circumvent the problem are outlined below.
2.5 (a) Full-field initialization
A diagnosis and discussion of the problems of adiabatic initialization has been given by Wergen (1987), who also described the diabatic non-linear normal-mode initialization scheme introduced operationally at ECMWF in 1982. This scheme was used in conjunction with the Centre's Optimum Interpolation analysis until the latter was replaced by a variational data assimilation scheme early in 1996. The diabatic initialization scheme was modified in 1986 to include the treatment of tides specified below.
The approach involves estimating a steady large-scale diabatic forcing. This is applied as a tendency term that does not change during the iterations of the non-linear normal-mode initialization, thereby avoiding the problem of non-convergence. The tidal problem is dealt with by subtracting an estimate of the tidal tendency from the adiabatic tendency. The steps involved are:
| • Performing a short forecast (~2h) from the uninitialized analysis, using all parametrizations but suppressing the diurnal cycle in the radiation scheme; |
| • Time-averaging the diabatic tendencies computed at each timestep of this forecast; |
| • Projecting the time-averaged tendencies onto gravity modes, keeping only low-frequency components (with periods>11h); |
| • Subtracting tidal tendencies computed from a (10-day) time series of the most recent analyses; |
| • Adding the resulting fixed filtered tendencies to the adiabatic tendencies computed for each iteration of the initialization procedure; |
| • Computing the initialization increment in the usual way, using the modified tendencies from the preceding step. |
This approach, though rather cumbersome, worked well in practice.
2.5 (b) Incremental initialization
The simpler alternative approach of incremental initialization (Puri et al., 1982; Ballish et al., 1992) was adopted when ECMWF moved to the variational data assimilation scheme. This approach does not necessarily rely on using the normal-mode technique for initialization. In particular, it can be used in conjunction with the digital filtering technique described later.
We illustrate the approach here in its basic form. Let
denote a background state, and as before 
the uninitialized state and 
the initialized state. Let 
denote the result of an adiabatic initialization of 
. The incremental initialization scheme is then:
|
(64) |
Provided the differences between
and 
are small, the scheme preserves the large-scale
diabatic balance present in the background field. However, if a component
is handled poorly by the initialization procedure and is improved in 
compared with 
, the thermal tide for example, the improvement in 
may not carry through to the initialized analysis 
.As applied in a conventional data assimilation system, the state
is the short-range forecast that provides the background
for the analysis, 
is the result of the analysis and 
is the starting point for the next forecast. The approach has, however,
wider applicability. (64) could also be used, for example, to
provide initial conditions for a forecast from an analysis produced using
a different data assimilation system, such as that of a different forecasting
centre. In this case, 
would be the (uninitialized) analysis produced by a standard data
assimilation using the forecast model for which initialization is required,
and 
would be the analysis produced by the
different data assimilation system.Training Course Notes Front Page >>
Table of contents >>
Next Section >>
Previous Section >>
 |
07.06.2002 |
 |
© ECMWF |
