Initiation of ray tracing models: evolution of small-amplitude gravity wave packets in non-uniform background

Abstract

This study introduces a new computational scheme for the linear evolution of internal gravity wave packets passing over strongly non-uniform stratifications and background flows as found, e.g., near the tropopause. Focusing on linear dispersion, which is dominant at small wave amplitudes, the scheme describes general wave superpositions arising from wave reflections near strong variations of the background stratification. Thus, it complements WKB theory, which is restricted to nearly monochromatic waves but covers weakly nonlinear effects in turn. One envisaged application of the method is the formulation of bottom-of-the-stratosphere starting conditions for ray tracing parameterizations that follow nonlinear gravity wave packets into the upper atmosphere. A key feature in this context is the method’s separation of wave packets into up- and downward-propagating components. The paper first summarizes a multilayer method for the numerical solution of the Taylor–Goldstein equation. Borrowing ideas from Eliassen and Palm (Geophys Publ 22:1–23, 1961), the scheme is based on partitioning the atmosphere into several uniformly stratified layers. This allows for analytical plane wave solutions in each layer, which are matched carefully to obtain continuously differentiable global eigenmode solutions. This scheme enables rapid evaluations of reflection and transmission coefficients for internal waves impinging on the tropopause from below as functions of frequency and horizontal wavenumber. The study then deals with a spectral method for propagating wave packets passing over non-uniform backgrounds. Such non-stationary solutions are approximated by superposition of Taylor–Goldstein eigenmodes. Particular attention is paid to an algorithm that translates wave packet initial data in the form of modulated sinusoidal signals into amplitude distributions for the system’s eigenmodes. With this initialization in place, the state of the perturbations at any given subsequent time is obtained by a single superposition of suitably phase-shifted eigenmodes, i.e., without any time-stepping iterations. Comparisons of solutions for wave packet evolution with those obtained from a nonlinear atmospheric flow solver reveal that apparently nonlinear effects can be the result of subtle linear wave packet dispersion.

Appendices

Convergence proof for the multilayer method with the $C^1$-approach

This section shows a proof that the multilayer method using the matching conditions from the $C^1$-approach (54) converges to the solution of Eq. (51) for fixed $\omega$. It relies on writing the method as a one-step method for the equivalent first-order system

$$\begin{aligned} \begin{pmatrix} w\\ q \end{pmatrix}^\prime = \begin{pmatrix} 0 &{}\quad 1 \\ -m^2(z) &{}\quad 0 \end{pmatrix} \begin{pmatrix} w\\ q \end{pmatrix}. \end{aligned}$$

(62)

In particular, given a point $z_0$, initial conditions $w_0, q_0$ and a step size $\Delta _z$, we solve the constant-coefficient system

$$\begin{aligned} \begin{pmatrix} w^{(1)}\\ q^{(1)} \end{pmatrix}^\prime = \begin{pmatrix} 0 &{}\quad 1 \\ -m^2(z_0 +\frac{\Delta _z}{2}) &{} \quad 0 \end{pmatrix} \begin{pmatrix} w^{(1)}\\ q^{(1)} \end{pmatrix}, \end{aligned}$$

(63)

due to the initial conditions $w^{(1)}(z_0)=w_0, q^{(1)}(z_0)=q_0)$ and set $w_1 = w^{(1)}(z_0 + \Delta _z), q_1 = q^{(1)}(z_0 + \Delta _z)$, which are the initial conditions for the next step. Since $q=w^\prime$, owing to the first row of the equation system, the global solution and its first derivative are continuous. Lara [18] shows that such one-step methods for first-order systems converge to the solution of the corresponding system. But the one-step method from Eq. (63) is equivalent to the multilayer method, as the eigenvalues of the matrix are $\pm i m(z_0 + \frac{\Delta _z}{2})$ and the solution is a linear combination of the exponentials of the eigenvalues, which is exactly the same as for the multilayer method. This puts the use of the $C^1$-approach on a solid foundation.

Grid convergence study for the multilayer method for wave packets

This section is dedicated to a grid convergence study for the multi layer method. We chose two test cases already presented in the main part of this paper, namely the propagation of a wave packet through a realistic tropopause profile, as shown in Fig. 5, and a wave packet refracted by a weak background wind, i.e., case REFR from Sect. 4 with the same setups for initial wave packets as in the respective test cases.

The reference solution $w_\mathrm{ref}$ is, in both cases, computed with $2^9$ layers, and we compare it to the solution $w_n$ computed with $n = 2^2, 2^3, \ldots , 2^8$ layers. The time point of the comparison was chosen such that the initial wave packet has passed the area of non-uniform background. The real solution was computed with a vertical resolution of $50\,\mathrm {m}$. We compute the absolute error between those 2 solutions in the $l_\infty$-norm, i.e.,

$$\begin{aligned} l_\infty (n) = \max \left| w_n -w_\mathrm{ref} \right| , \end{aligned}$$

(64)

but we normalize this error by the maximum vertical wind speed of the reference solution $\max \left| w_\mathrm{ref} \right|$.

The results are displayed in Fig. 14. For the wave packet moving through a tropopause, which can be seen in the left panel, we observe that the error decreases with a fixed convergence order. The slope of a least square fit is $-2.07$, which indicates that we have a quadratic convergence.

For the wave packet in the REFR case, we see that a certain threshold in layers has to be surpassed (which seems to be somewhere between 32 and 64 layers) in order to obtain acceptable results. But once this requirement is met, we can also see that the error decreases. A least squares fit here has a slope of $-1.52$, and when considering only the data points from 32 layers onward, we have a slope of $-2.23$, which lets us conclude that the convergence order is 2.

Fig. 14

Error development for increasing number of layers for a wave packet moving through the tropopause profile (48) (left) and for the case REFR (right)