## Abstract

Gravity waves arise in gravitationally stratified compressible flows at low Mach and Froude numbers, and these waves impose a sharp restriction on the time step. This paper presents a filtering strategy for the fully compressible equations based on normal-mode analysis that is used throughout the simulation to compute the fast dynamics and is able to damp only chosen modes. This method is based on an asymptotic analysis and respects the dynamics of gravity waves for thin layers. Finally, the filtering method is tested on a series of examples.

## 1. Introduction

Gravity waves are present in gravitationally stratified compressible flows at low Mach and Froude numbers, such as atmospheric flows and stars. These waves present a particular challenge when considering an adaptive mesh refinement (AMR) framework. Most of the models developed for single-grid calculations are hydrostatic, where the pressure evolution equation is replaced by the hydrostatic equation. In an AMR calculation, where the computational domain contains several grid sizes and the time step on the finer grids is subcycled (e.g. Jablonowski 2004), the formulation must be non-hydrostatic. The presence of the pressure evolution equation allows one to match the pressures easily at the coarse–fine interfaces in a time subcycled simulation and to avoid non-physical effects. However, the time step is reduced because the Courant–Friedrichs–Lewy (CFL) condition now depends on the velocity of the fastest gravity wave or even the fastest acoustic wave if the acoustic dynamics is not filtered. Therefore, an important idea in devising methods for atmospheric flows in an AMR framework is to identify each of the time scales and to provide a means of separating them. In Gatti-Bono & Colella (2006), we presented an all-speed projection method that allowed the elimination of the CFL restriction related to the acoustic waves.

In this paper, we summarize an all-speed projection method that allows the advection dynamics to be separated from the acoustic dynamics. Then, we propose a method to isolate each of the gravity wave modes that are faster than the underlying fluid advection, and therefore propose a method where the time step is constrained only by the advective motions.

## 2. An all-speed formulation

We consider an inviscid compressible fluid described by the compressible Euler equations in perturbational form:
2.1
2.2
and
2.3
Here, *ρ*_{o}(*z*) is the hydrostatic density, is the perturbational density, *p*_{o}(*z*) is the hydrostatic pressure defined as
2.4
and is the perturbational pressure.

These equations contain different time scales, i.e. acoustic waves, gravity waves and advection, which are intertwined in the current form of the Euler equations. We introduce a splitting of velocity and pressure as a first step to separate out the different time scales:
2.5
and
2.6
We first give the names of the different quantities from the splitting and then give their definitions in the next paragraph: **u** is the total velocity, **u**_{d} is the anelastic velocity, **u**_{p} is the curl-free velocity, **u**_{h} is the harmonic velocity, *π*_{H} is the ‘hydrostatic-like’ pressure associated with the gravity waves, *π*_{I} is the pressure associated with the incompressible effects, *ψ* is the pressure that corresponds to the vertical motions generated by the horizontal gravity waves and *δ* is the acoustic pressure.

The velocities **u**_{d} and **u**_{p} can be obtained through projections of the total velocity and they are defined as
2.7
and
2.8
where
2.9
2.10
and
2.11
They satisfy the relationships
2.12
2.13
and
2.14
The velocity **u**_{h} is a space-dependent quantity, but it is not time dependent and is computed once for the initial conditions.

We have introduced a parameter *η*_{o}:
2.15
The choice of *η*_{o} is detailed in Gatti-Bono & Colella (2006), and it is driven by the requirement of the splitting that the gravity wave dynamics separates from the acoustic dynamics. We note that *η*_{o} corresponds to the isentropic density associated with the pressure distribution *p*_{o}. With this definition of *η*_{o}, equation (2.12) is a constraint similar to the pseudo-incompressible constraint introduced by Durran (1989). However, in the current work, the constraint is imposed only on the anelastic part of the velocity, whereas it was imposed on the total velocity in Durran (1989).

A decomposition of the pressure must be done along with the splitting of the velocities. The different pressures in equation (2.6) are defined below: 2.16 2.17 2.18 and 2.19 with 2.20 and 2.21

Using the splitting, equations (2.1)–(2.4) become
2.22
2.23
2.24
and
2.25
where *N* is the Brunt–Väisälä frequency defined by
2.26

In this form, equations (2.22) and (2.23) contain the anelastic dynamics and equations (2.24) and (2.25) contain the acoustic dynamics. The next step is to isolate the normal modes.

## 3. Isolating the normal modes

Equations (2.22) and (2.23) can be cast into a system of hyperbolic equations: 3.1 and 3.2 where 3.3 with 3.4 3.5 and 3.6

Gatti-Bono & Colella (2006) showed that, by using an asymptotic analysis along with a change of variables, the system of equations (3.1) and (3.2) becomes
3.7
and
3.8
where *ϵ* is the aspect ratio of the domain and with
3.9
3.10
and
3.11

When *χ* is negative, the operator is positive definite and equations (3.1) and (3.2) describe the motion of an infinite collection of horizontally propagating gravity waves, one for each eigenmode of the second-order self-adjoint operator . The operator has eigenvectors *r*^{k} and eigenvalues *λ*^{k}, and the eigenvalues *λ*^{k} are related to the gravity wave speeds *c*^{k} by
3.12

Figure 1 shows the gravity wave speeds for the first example presented in §6*a*(ii) along with the first four eigenmodes. As expected, only a few modes travel with a speed larger than a typical mean fluid velocity of 20 m s^{−1} (solid line) and the number of sign changes in the eigenmode increases as the mode number increases.

When *χ* is positive, the first few eigenvalues can become negative, and this is a direct consequence of including the potential dynamics into the asymptotic analysis. We observe that this has a stabilizing effect, as the modes with negative eigenvalues do not support gravity waves. A detailed discussion is available in Gatti-Bono & Colella (2006).

If the domain is not regular, we introduce a change of coordinates that maps the coordinates *z* of the computational domain with embedded boundaries into coordinates *z** on a regular grid:
3.13
where *L* is the height of the domain and *z*_{o}(*x*) is the mountain profile. The operator becomes
3.14
where *η*_{o}, *ρ*_{o}, *ξ*, *χ* and *ζ* are taken in the new reference frame. The operator is still positive definite when *χ* is negative, and we can therefore compute eigenvectors and eigenvalues for this operator.

## 4. An all-speed algorithm with filtering

We discretize in time using Δ*t*, and we solve for the total velocity **u**=(*u*,*w*), the anelastic velocity **u**_{d}=(*u*_{d},*w*_{d}), the curl-free velocity **u**_{p}=(*u*_{p},*w*_{p}), the total pressure *p*, the acoustic pressure *δ*, the auxiliary pressures *ψ*, *π*_{H} and *π*_{I}, the density *ρ* and the perturbational density .

First, we solve for and : 4.1 and 4.2 where 4.3 Note that and depend on only the variables

*x*and*t*. The * represents variables in the mapped coordinate system when there is topography.Then, the anelastic velocity is partially advanced in time using only the advection terms in equation (2.23): 4.4 We compute the quantities at half time steps with a Godunov method to form

**A**_{d}**u**^{n+1/2}.We use equation (2.22) to temporarily advance the density variables

*ρ*and : 4.5 4.6 and 4.7Equation (2.16) is used to temporarily advance

*π*_{H}: 4.8 Here,*L*is the vertical upper bound of the domain.We filter out fast modes in

*π*_{H}in the mapped coordinate space without embedded boundaries: 4.9 with given in §5*a*(ii). From : we compute : 4.10We filter out the fast modes in : 4.11 and 4.12

Using equation (2.17), we solve for : 4.13

Using equation (2.18), we solve for

*ψ*^{n+1}: 4.14Combining equations (2.24) and (2.25), we implicitly advance the acoustic pressure: 4.15 where 4.16

Then, we advance the curl-free velocity using equation (2.24): 4.17

The perturbational density is advanced using equation (2.22): 4.18

Then

*π*_{H}is updated: 4.19We filter out fast modes in

*π*_{H}in the mapped coordinate space without embedded boundaries: 4.20 From , we compute .We filter out the fast modes in : 4.21 and 4.22

Then

*ψ*is updated: 4.23We finish advancing the divergence-free velocity: 4.24

We compute the fully filtered divergence-free velocity: 4.25 with given in §5

*b*. We then compute in the space of embedded boundaries.We project the divergence-free velocity: 4.26

## 5. Details of the filtering implementation

In the algorithm of the previous section, we need to filter *π*_{H} and *u*_{d}. This section details the equations and boundary conditions used for the filtering.

### (a) Projections and filtering for *π*_{H}

We solve two systems of equations, each with their own set of boundary conditions.

#### (i) System with homogeneous boundary conditions

For this system,
5.1
and
5.2
with boundary conditions
5.3
and
5.4
This problem comes down to finding the eigenvectors *r*^{k} and eigenvalues of with boundary conditions
5.5
and
5.6

#### (ii) System with inhomogeneous boundary conditions

For this system, 5.7 with boundary conditions 5.8 and 5.9 The total pressure can be written as 5.10 and 5.11

### (b) Projections and filtering for *u*_{d}

We write *u*_{d} as
5.12
where
5.13
with
5.14
and
5.15

## 6. Results

### (a) Perturbational gravity wave test problem

In Gatti-Bono & Colella (2006), we set up an example to verify numerically that the algorithm and the splitting of the full compressible equations allow us to recover results from the asymptotic analysis. We initialized the different variables so that, initially, we have a Gaussian distribution in selected modes only. We showed that, if the initial data are only in the fastest mode, the wave propagated only in that mode when we used a time step that was smaller than the one dictated by the CFL condition associated with the fastest gravity wave. Here, we add the filtering algorithm as discussed above, and we show that the computation is stable for time steps that are larger than the one dictated by the CFL condition associated with the fastest gravity wave. We also use examples in which the initial data are not in the fastest (filtered) mode.

#### (i) Summary of the initialization

The solenoidal velocity is initialized as
6.1
where is a selection of mode numbers, *r*^{k} is the eigenvector in the *k*th mode, *G*_{k}(*x*) is a Gaussian function
6.2
and the pressure *π*_{H} is
6.3
where *c*^{k} is the speed of the *k*th gravity wave.

The remaining variables are initialized according to their definitions or the equations of motion (see Gatti-Bono & Colella 2006).

#### (ii) Example 1: initial data in the fastest mode only

In this example, the initial data are chosen so that they are only in the fastest mode: *k*_{1}=0 and *p*=1. We use different domain sizes that correspond to different aspect ratios *ϵ* and different grid refinements, as shown in table 1.

Figure 2 shows the projection onto the fast mode of the horizontal velocity and of the pressure *π*_{H} for different time steps with and without filtering. The fastest (respectively, second fastest) gravity wave propagates with a speed of 164 m s^{−1} (respectively, 54.7 m s^{−1}), and, to match the CFL conditions for an aspect ratio of *ϵ*=0.0125 and a grid of 128×64, we have to use a time step smaller than 48.7 s (respectively, 146 s). Figure 2 shows that the filtering preserves the shape and propagation speed, but introduces some damping in the solution, which is expected because a Helmholtz solve is used to compute the one-dimensional fast solution. However, it should be noted that the filtering of just this one mode alone allows us to multiply the time step by a factor of 3 without losing important features of the solution, whereas the unfiltered solution becomes unstable after just a couple of time steps. Also note that this is an artificial example that focuses on the fastest mode of the solution when, for most mesoscale atmospheric modelling, the few fastest gravity waves may not be significant for the dynamics but may restrict the time step significantly. Therefore, ultimately, we do not expect that some damping in the fastest mode would generate accuracy issues for mesoscale simulations.

Figure 3 shows that the solution is in the asymptotic regime for an aspect ratio of *ϵ*=0.0125 and that, for this aspect ratio, the damping diminishes when the grid is refined.

#### (iii) Example 2: initial data in modes 2 and 3

Here, we use two different initial data. For case (i), the initial data are only in the second mode (i.e. *k*_{1}=2 and *p*=1); and for case (ii), they are a combination in the second and third modes (i.e. *k*_{1}=2, *k*_{2}=3 and *p*=2). Figure 4 shows the projections of the divergence-free velocity and pressure *π*_{H} onto the second mode for case (i) (top) and for case (ii) (middle) and onto the third mode for case (ii) (bottom). These graphs show that the filtering method has only a negligible effect on the accuracy of the solution at a fixed time step when the filtered mode is not present in the initial data. Unlike the results in the previous section, filtering does not introduce damping. When the time step is multiplied by a factor of 3, very little damping is introduced and the solution is almost identical to that for a smaller time step that satisfies the CFL condition for the fastest gravity wave.

### (b) Mountain lee waves

This is the classic example of mountain lee waves that can be found in Durran (1986). A detailed description of the initial conditions and parameters can be found in Gatti-Bono & Colella (2006). In this example, a uniform wind flows from left to right at a speed of 20 m s^{−1} over a 600 m high mountain. The mountain is described as a ‘witch of Agnesi’,
6.4
with a width *a* of 10 km, a height *h* of 600 m and a centre *x*_{o} positioned at 72 km. The atmosphere is a two-layered atmosphere with a constant Brunt–Väisälä frequency *N* of 0.01 or 0.02 s^{−1}.

Figure 5 shows the mountain lee waves example run with and without filtering at different time steps for the unstable cases, i.e. the cases where gravity waves are expected. The filtering allows us to double the time step but damps the gravity wave response because this response is essentially in the first mode. However, the main features of the solution are retained when filtering is used. It is interesting to note that this example shows that the filtering method can be used successfully when an embedded boundary is present and in a real atmospheric setting with a significant gravity wave response present.

## 7. Conclusion

This paper presents a method to filter gravity waves out of the fully compressible equations, which is based on the normal-mode analysis. The scheme was embedded in an all-speed algorithm. A series of examples were presented to test the filtering method on the algorithm presented in Gatti-Bono & Colella (2006). The method was shown to damp only the modes that are being filtered, keeping the main features of the solution, such as wavelengths, unaltered. It was also run successfully on the mountain lee wave example with embedded boundaries. In both sets of examples, the computation with filtering was stable when the time step was multiplied by a factor of 2 or 3, whereas the computation without filtering would become unstable after only a few time steps.

## Acknowledgements

We would like to acknowledge Dr Joseph Tribbia for suggesting that we look at the vertical normal-mode analysis.

## Footnotes

One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.

- © 2009 The Royal Society